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Abstract 

A  Non-Linear  Finite  Element  Model  for  the  Determination  of  Elastic 
and  Thermal  Properties  of  Nanocomposites 

by 

Paul  Elsbernd 


This  thesis  presents  a  non-linear  model  for  the  thermal  and  elastic  properties  of  single- 
walled  carbon  nanotube  reinforced  polymer  composites.  Finite  Element  Analysis  (FEA), 
in  conjunction  with  the  Embedded  Fiber  Method  (EFM),  is  used  to  calculate  the  effective 
stress-strain  curve  and  thermal  conductivity  of  the  composite  material.  First,  the 
geometry  of  a  user-defined  volume  fraction  of  nanotubes  is  randomly  generated  and  their 
properties  are  incorporated  into  the  polymer  matrix  using  the  EFM.  Non-linear  FEA  is 
next  performed  to  account  for  the  non-linear  properties  of  the  polymer  matrix  and  the 
carbon  nanotubes.  Finally,  Monte  Carlo  Analysis  of  five  hundred  random 
micro  structures  is  performed  to  capture  the  stochastic  nature  of  the  random  fiber 
generation  and  to  derive  statistically  sound  results.  The  model  is  validated  by 
comparison  with  several  different  experiments  reported  in  the  open  literature. 
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Chapter  1 

Introduction 


1.1  Motivation 

As  nanotechnology  has  blossomed  over  the  last  two  decades,  countless  applications  for  a 
variety  of  nanoparticles  have  been  demonstrated  or  proposed.  The  use  of  gold 
nanoparticles  to  destroy  tumors  [1],  of  carbon  nanotubes  to  improve  bullet  proof  vests 
[2],  and  of  fullerene  nanoparticles  to  detect  biological  markers  of  disease  [3]  are  a  few 
examples  that  encompass  only  a  tiny  fraction  of  the  countless  possible  applications  of 
nanoparticles.  Arguably,  the  most  commonly  proposed  application  of  nanoparticles  is 
their  use  in  composite  materials.  While  the  choice  of  particle  strongly  depends  on  the 
functionality  sought  for  the  composite,  one  of  the  most  widely  studied  nanocomposites  is 
the  carbon  nanotube  reinforced  composite. 

Carbon  nanotubes  offer  exceptional  mechanical,  thermal,  and  electrical  properties  at  an 
almost  insignificant  weight  on  the  macroscale.  Theoretically,  single-walled  carbon 
nanotubes  (SWCNTs)  exhibit  elastic  moduli  of  the  order  of  1  TPa  and  fracture  strains  of 
10-30%  [4-8],  Initial  experiments  did  not  identify  such  a  specific  value  for  the  elastic 
modulus  due  to  undeveloped  experimental  techniques  and  carbon  nanotubes  that  were  not 
perfectly  prepared  or  treated.  Troiani  et  al.  [9]  reported  an  elastic  modulus  of 
305  +  lSGGPa.  On  the  other  extreme,  Treacy  et  al.  [10]  report  1JS +0.9  TPa.  Recent 
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improvements  in  nanotube  fabrication  and  new  experimental  techniques  have  yielded 
elastic  properties  that  are  more  consistent  with  theoretical  predictions.  Wong  et  al.  [11] 
and  Yu  et  al.  [8]  reported  elastic  moduli  of  1.28+  0.59  TPa  and  0.95  TPa,  respectively. 
These  findings  indicate  that  SWCNT’s  have  elastic  moduli  that  are  three  times  that  of 
carbon  fibers,  and  five  times  that  of  steel  at  one  sixth  the  weight  [12], 

In  addition  to  their  outstanding  mechanical  properties,  carbon  nanotubes  exhibit 
exceptional  thermal  and  electrical  properties.  The  theoretical  thermal  conductivity  of 

fir* 

SWCNTs  is  usually  reported  as  6,000  —  at  room  temperature  [13],  which  is 

approximately  three  times  the  thermal  conductivity  of  diamond,  and  approximately  10* 
times  that  of  most  polymers  [14],  This  value  is  in  fact  highly  temperature  dependent  as 
reported  by  Grujicic  et  al.  [15],  ranging  from  values  as  high  as  12,000  —  at  100K  to 

771— K 

W 

about  2,000  at  300K.  Carbon  nanotubes  are  also  found  to  have  exceptional 
electrical  conductivity,  of  the  order  of  104  ^107  — ,  approximately  20  decades  higher  than 
that  of  most  polymers  [14,16]. 

These  outstanding  properties  make  carbon  nanotubes  extremely  interesting  for  a  variety 
of  applications.  Individually,  applications  for  CNTs  include  field  emission  sources, 
Lithium-ion  batteries,  capacitors,  actuators,  molecular  sensors,  gas  and  hydrogen  sources, 
scanning  probe  tips,  and  countless  electronics  applications  [17,18],  However,  of  more 
interest  to  this  current  study  are  the  potential  applications  for  nanotubes  in  composite 
materials;  polymer  composites  in  particular.  Their  mechanical  properties  make  them 
ideal  candidates  to  use  as  reinforcing  agents,  either  alone  or  as  filler  with  other 
reinforcements,  such  as  carbon  or  glass  fibers.  Their  thermal  properties  make  them 
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desirable  for  thermal  management  applications,  and  their  electrical  properties  make  them 
ideal  for  numerous  other  composite  applications  such  as  electromagnetic  shielding, 
preventing  electrostatic  charging,  and  damage  sensing  [19,20].  Collectively,  these 
properties  make  nanotube  composites  very  interesting  to  many  different  fields  including 
the  aerospace,  automotive,  medical,  and  defense  industrial  sectors. 

Because  of  the  significant  promise  in  these  new  composites,  a  need  to  predict  their 
physical  properties  has  developed.  It  is  in  this  context  that  the  motivation  for  the  current 
work  is  derived.  The  capacity  to  accurately  predict  properties  of  nanocomposites  is 
particularly  needed  because  the  sheer  number  of  different  possible  composites  limits 
experimental  characterization.  Not  only  are  there  countless  different  polymers  to  be  used 
as  the  matrix  material  of  the  composite,  but  the  effect  of  various  volume  fractions  of 
nanotubes  also  needs  to  be  studied  for  each  composite.  These  two  factors  are  the  primary 
motivations  for  an  accurate  and  versatile  model  of  carbon  nanotube  composites  behavior. 

1.2  Obstacles  to  Accurate  Experimentation 


Besides  the  overwhelming  number  of  possible  composites,  there  are  four  major  issues 
influencing  the  experimental  accumulation  of  accurate  information  about  nanotube 
composites.  They  are:  inadequate  interfacial  bonding,  poor  dispersion,  waviness  and 
alignment  issues,  and  nanoscale  measurement  concerns.  Interfacial  bonding  describes  the 
ability  of  the  carbon  nanotube  to  chemically  bond  to  the  matrix  material.  This  point 
relates  to  the  inherent  physical  structure  of  carbon  nanotubes.  Pristine  nanotubes  are 
essentially  defect-free  and  atomically  smooth.  This  fact  inhibits  interfacial  bonding  with 
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the  matrix  structures  in  which  CNTs  are  mixed  and  thus  tremendously  inhibits  their 
reinforcing  capabilities. 

This  problem  has  been  approached  by  several  researchers  and  the  most  effective  solution 
has  been  functionalization.  Surface  functionalization  is  the  process  by  which  various 
chemical  functional  groups  are  attached  to  the  nanotube  sidewalls  and/or  endcaps  using 
either  covalent  or  non-covalent  bonding.  This  process  was  first  proposed  for  nanotubes 
by  Tiano  et  al.  [21].  By  specifically  choosing  functional  groups  based  upon  the 
composite  matrix  material,  a  dramatic  increase  in  the  interfacial  bonding  in  the  composite 
can  be  achieved.  Zhu  et  al.  [22]  report  that  by  utilizing  acid  treatment  and  fluorination  of 
the  nanotubes,  they  were  able  to  achieve  a  30%  increase  in  the  elastic  modulus  of  an 
epoxy  composite  with  the  addition  of  only  1%  weight  fraction  (wt.)  of  SWCNTs. 
Functionalization  is  an  excellent  method  because  it  allows  one  to  modify  nanotubes  with 
functional  groups  that  offer  the  best  bonding  for  any  specific  matrix  material,  and  it  only 
has  a  small  effect  on  the  properties  of  the  nanotube  itself.  Garg  and  Sinnot  [23]  report 
that  even  a  high  degree  of  sidewall  functionalization  will  degrade  the  mechanical  strength 
of  SWCNTs  by  only  15%. 

The  second  obstacle  is  that  of  dispersion.  Due  to  the  extraordinarily  large  aspect  ratio  of 
SWCNTs,  up  to  1,000,  nanotubes  tend  to  agglomerate  to  form  ropes  or  bundles  by  means 
of  intrinsic  Van  der  Waals  attraction  [24],  One  rope  can  contain  hundreds  of  nanotubes, 
but  the  nanotubes  can  readily  slide  relative  to  each  other  due  to  the  low  shear  modulus  of 
the  bundles  [25-27],  This  fact  prevents  good  dispersion  of  nanotubes  in  the  matrix  and 
not  only  limits  the  reinforcing  capabilities  of  the  CNTs,  but  it  can  actually  degrade  the 
initial  properties  of  the  bulk  matrix  material  alone.  Many  methods  have  been  studied  to 
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address  this  problem  including  ultrasonication,  high  shear  mixing,  surfactants, 
functionalization,  and  several  others  [22],  The  most  successful  method  seems  to  be  a 
combination  of  the  above  methods  with  functionalization  playing  a  primary  role. 

In  addition  to  dispersion  and  bonding  issues,  waviness  and  alignment  are  other  major 
challenges  faced  by  experimentalists.  A  large  number  of  studies  have  been  published 
researching  methods  of  better  aligning  nanotubes  in  a  composite.  The  methods  of 
electrospinning,  magnetic  field  inducing,  liquid  crystal  inducing,  shear  flow,  extrusion  or 
ejecting,  melt  spinning,  and  mechanical  stretching  are  among  many  that  have  been 
explored  [17,22].  Furthermore,  lone  carbon  nanotubes,  even  in  composites,  are  generally 
assumed  to  be  straight  and  are  often  modeled  as  such.  In  actuality,  only  quite  short 
nanotubes,  less  than  about  100  ffm,  are  found  to  be  straight  or  “quasi-straight”  in 
composites  [20],  In  fact,  when  normal  nanotubes  are  incorporated  into  composites,  the 
bonding  with  the  matrix  material  in  conjunction  with  curing  of  the  composite,  in  the  case 
of  epoxies,  causes  waviness  to  be  observed.  This  is  due  to  the  low  bending  rigidity  of 
carbon  nanotubes,  of  the  order  of  1  X  10” 25  —  [28].  Additionally,  it  has  been  shown 

that  current  nanotube  fabrication  methods  result  in  nanotubes  with  a  large  distribution  of 
diameters  and  lengths  [29,30].  This  fact  will  in  turn  lead  to  a  distribution  of  waviness 
exhibited  in  the  nanotubes,  increasing  with  nanotube  length.  Most  current  models 
assume  that  the  nanotubes  in  the  composite  are  straight  due  to  the  computational  and 
theoretical  difficulties  that  arise  in  accounting  for  this  parameter  in  their  approach.  Only 
a  few  authors,  including  Shi  et  al.  for  mechanical  properties  [31],  Li  et  al.  for  electrical 
properties  [19],  and  Zhang  et  al.  for  thermal  properties  [32],  have  incorporated  waviness 
into  their  models.  Neither  of  these  authors  studied  the  mechanical  properties  of 
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composites.  The  model  proposed  in  this  paper  includes  both  waviness  and  alignment  in 
its  calculation  of  thermal  and  mechanical  properties  of  the  composite.  This  is  one  of  the 
primary  advantages  of  the  proposed  model  over  other  current  models. 

The  final  obstacle  relating  to  the  accumulation  of  knowledge  about  the  effective 
properties  for  different  composites  is  the  scale  on  which  experiments  must  be  conducted. 
Though  perhaps  less  challenging  than  observing  nanotube  properties  directly, 
nanocomposite  experimentation  is  difficult  because  it  is  both  costly  and  time  intensive 
due  to  the  equipment  and  methods  that  must  be  used.  While  nanotubes  are  costly  by 
themselves,  the  fabrication  of  the  composites  is  quite  difficult  and  time  consuming  as 
well.  For  polymer  composites,  which  are  of  considerable  interest,  most  researchers  first 
functionalize  the  nanotubes,  disperse  them  in  an  epoxy  matrix  material,  and  then  cure  the 
resin.  Zhu  et  al.  reported  this  process  to  take  a  minimum  of  18  hours  [22],  Any  misstep 
in  the  process  would  require  starting  over.  Then,  to  perform  an  analysis  on  the  sample, 
ATR-FTIR,  SEM/EDAX,  optical  microscopes,  DMA  systems,  and  other  systems  are 
used  [22,33],  Finally,  mechanical  testing  is  performed  in  a  variety  of  ways.  The 
complexity,  expense,  and  time  required  for  this  testing  are  some  of  the  primary 
motivations  for  a  general  model  of  nanotube  composites  to  be  developed. 


1.3  Current  CNT  Composite  Models 


In  pursuing  a  versatile  model,  many  different  models  have  been  developed  with  a  variety 
of  unique  approaches.  Tserpes  et  al.  [34]  developed  a  multi-scaled  modeling  approach  to 
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test  the  tensile  behavior  of  composites  in  which  nanotubes  are  modeled  independently, 
then  approximated  as  a  beam  element,  and  integrated  into  a  representative  volume 
element  (RVE)  which  can  eventually  be  solved  using  finite  element  analysis  (FEA).  Li 
and  Chou  [35]  presented  another  multi-scale  approach  to  test  compressive  behavior. 
Xiao  and  Gillespie  Jr.  [36]  reported  a  nanomechanics  model  for  predicting  elastic 
properties  of  SWCNTs  as  composite  reinforcement,  and  then  used  a  micro-mechanics 
model  to  calculate  the  elastic  properties  of  the  composite  as  a  whole.  Shi  et  al.  [31]  used 
a  micromechanics  model  as  well,  and  include  waviness  and  agglomeration  in  their 
calculation  of  the  elastic  properties.  Song  and  Youn  [37]  calculated  the  effective  thermal 
conductivity  of  nanotube/polymer  composites  using  the  control  volume  finite  element 
method.  All  of  these  models  use  a  single  carbon  nanotube  that  traverses  the  entire  RVE. 
They  assume  that  the  effective  properties  of  the  composite  as  a  whole  can  be  calculated 
based  on  a  representative  RVE  with  a  single  nanotube,  primarily  due  to  computational 
simplicity.  This  assumption  is  contended  in  this  paper  and  is  not  made  by  the  current 
model. 

Beyond  the  models  that  calculate  effective  properties  based  upon  a  single  nanotube  RVE, 
there  are  several  models  that  do  not  make  this  assumption.  Xue  [38]  presented  a 
numerical  method  based  on  Maxwell  theory  to  calculate  the  thermal  conductivity  of 
nanotube/oil  and  nanotube/decene  composites.  Bagchi  and  Nomura  [39]  approximated 
nanotubes  as  spheroidal  inclusions  and  used  effective  medium  theory  to  calculate  the 
thermal  conductivity  of  aligned  nanotube/polymer  composites.  One  of  the  most  versatile 
models  is  presented  by  Zhang  and  Tanaka  [32],  Using  the  hybrid  boundary  node  method 
in  conjunction  with  the  Fast  Mulitpole  Method,  they  calculated  the  thermal  properties  of 
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CNT  composites  while  including  waviness  and  alignment  for  a  variety  of  nanotube 
configurations.  Odegard  et  al.  [40]  modeled  the  nanotubes  as  an  effective  continuum 
fiber  using  the  equivalent-continuum  modeling  method  along  with  a  micromechanics 
approach  to  calculate  mechanical  properties  of  CNT/polymer  composites.  Seidel  and 
Lagoudas  [41]  modeled  nanotubes  using  a  composite  cylinders  approach  and  used  the 
Mori-Tanaka  methods  to  obtain  the  effective  elastic  properties  for  the  composite  as  a 
whole.  They  considered  both  aligned  and  randomly  oriented  nanotubes  and  also  consider 
nano  tube  bundling.  Finally,  Li  and  Chou  presented  a  model  for  damage  sensing  by 
calculating  electrical  properties  of  CNT/glass-fiber  composites  using  the  finite  element 
method. 

Each  of  the  methods  detailed  above  provide  important  and  unique  contributions  to  the 
modeling  of  CNT-reinforced  composites,  but  there  is  one  major  aspect  that  nearly  every 
model  excludes;  non-linearity.  Tserpes  et  al.  [34]  are  the  only  known  authors  who 
reported  a  stress-strain  curve  for  their  composite  calculations,  but  they  considered  only  a 
single-nanotube  that  is  perfectly  aligned  and  traverses  the  entire  RVE.  Without  the 
consideration  of  the  non-linear  elastic  properties  of  the  polymer  matrix  and/or  those  of 
the  carbon  nanotubes  themselves,  a  stress-strain  relationship  cannot  be  obtained  for  the 
composite.  The  stress-strain  curve  is  one  of  the  most  important  properties  needed  to 
characterize  a  composite  and  is  thus  the  primary  motivation  and  result  of  the  model 
detailed  in  this  thesis. 

The  model  proposed  herein  reliably  simulates  CNT-reinforced  composite  properties  and 
provides  qualitative  deformation  information  regarding  the  composite  structure  as  a 
whole.  Esteva  [42]  reported  a  finite  element  model  in  which  nanotube-composite 
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structures  are  modeled  using  finite  elements  and  the  Embedded  Fiber  Method,  originally 
developed  to  model  steel  bar  reinforcements  in  concrete  [43-45].  This  useful  approach 
allows  a  simple  square  mesh  of  the  entire  RVE  to  be  generated,  circumventing  a 
complicated  meshing  process  while  still  producing  accurate  results.  In  ordinary  FEM 
routines,  human  intervention  must  take  place  to  locate  fiber-mesh  interfaces  and  to 
increase  the  mesh  refinement  in  these  areas.  The  simple  square  mesh  used  by  the  EFM 
allows  for  an  RVE  with  a  realistic  volume  fraction  of  nanotubes  and  the  addition  of  non¬ 
linearity  to  be  modeled  with  a  reasonable  amount  of  computational  effort.  This  model 
includes  the  consideration  of  nanotube  waviness  and  alignment  as  well.  This  concept, 
when  applied  to  nanotube-composites,  is  novel  and  effective.  Esteva  was  able  to 
calculate  the  effective  modulus  of  elasticity  and  thermal  conductivity  of  a 
nanotube/polymer  composite,  but  his  model  did  not  include  a  non-linear  stress-strain 
relationship  or  a  non-linear  thermal  conductivity.  The  proposed  model  builds  upon  the 
work  of  Esteva,  but  it  also  incorporates  the  non-linear  elastic  properties  of  polymer 
composites  using  an  incremental  approach,  and  studies  the  non-linear  thermal  properties 
of  SWCNTs  using  both  incremental  and  iterative  approaches.  The  non-linear  elastic 
approach  allows  the  calculation  of  stress-strain  relationships  along  with  descriptive 
deformation  behavior  before  and  after  the  onset  of  plastic  deformation  of  the  composite. 
The  non-linear  thermal  approach  allows  for  the  effect  of  non-linear  thermal  conductivity 
to  be  examined  and  its  relevance  in  modeling  to  be  discussed.  In  addition  to  the  non¬ 
linear  issues,  the  model  herein  proposes  the  introduction  of  several  realistic  parameters, 
including  length  and  diameter  distributions  of  SWCNTs  as  well  as  a  more  accurate  model 
of  nanotube  waviness.  Each  of  these  aspects  is  considered  to  make  the  model  versatile 
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with  the  ultimate  goal  of  producing  a  model  that  closely  agrees  with  experimental  results. 
In  the  numerical  results  section  of  this  thesis,  a  comparison  to  the  experimental  results  of 
Sun  et  al.  [46]  and  Zhu  et  al.  [17,22]  is  presented  for  the  stress-strain  relationships 
calculated  by  the  proposed  model.  The  effective  thermal  conductivities  calculated  by  the 
model  are  also  compared  to  six  different  experiments  [16,47-51]. 
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Chapter  2 

The  Representative  Volume  Element 

The  herein  proposed  approach  is  based  on  a  two-dimensional  finite  element  model  of  the 
composite  structure.  The  properties  of  the  composite  structure  as  a  whole  are  represented 
by  a  representative  volume  element  (RYE)  of  dimension  1  #tm  x  1  fiva  x  5  m.  The  small 
thickness  approximates  a  two  dimensional  structure  which  is  adopted  for  simplicity. 
There  are  three  distinct  steps  in  the  proposed  model;  geometry  generation,  fiber 
portioning,  and  finite  element  analysis.  First,  to  generate  the  nanotube  geometry  of  the 
RYE,  the  positions  and  orientations  of  nanotubes  in  the  RVE,  up  to  the  desired  volume 
fraction,  are  randomly  generated.  Second,  the  RVE  is  divided  into  a  square  mesh  of 
elements  and  the  fibers  are  partitioned,  or  divided,  according  to  which  element  they  fall 
in.  The  refinement  of  this  square  mesh  necessary  for  accurate  results  is  one  of  the  objects 
of  study  in  this  report.  Finally,  the  mechanical  or  thermal  properties  of  the  nanotubes  in 
each  element  are  integrated  into  the  composite  matrix  properties,  based  upon  which  kind 
of  study  is  being  performed,  and  non-linear  finite  element  analysis  is  utilized  to  find 
either  the  stress-strain  relationship  and  deformation  behavior,  or  the  thermal  conductivity 
of  the  embedded  fiber  composite. 
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2.1  Geometry  Generation 

The  first  step  in  accurately  modeling  any  composite  material  is  to  realistically  model  the 
geometry  of  the  reinforcing  agents.  Due  to  the  many  advances  that  have  been  made  in 
dispersing  nanotubes  and  breaking  up  ropes  and  agglomerates  [52,53],  as  well  as  the 
effects  that  these  factors  have  on  effective  composite  properties  [54],  a  geometry  of  lone 
nanotubes  is  assumed.  To  generate  an  accurate  geometry  within  the  RVE,  the  parameters 
of  nanotube  length,  diameter,  and  waviness  are  integrated  into  the  model. 

In  any  individual  batch  of  carbon  nanotubes,  varying  according  to  method  of  fabrication, 
lengths  of  individual  nanotubes  vary  anywhere  from  less  than  100  iwi  to  over  1  ffln. 
Both  Wang  et  al.  [29]  and  Ziegler  et  al.  [30]  reported  statistically  accurate  length 
measurements  of  a  large  number  of  nanotubes.  Wang  et  al.  reported  a  Weibull 
distribution  of  lengths  with  scale  and  shape  parameters  ofa=SM  IQ-*  and  b  =  2.4, 
respectively  [29],  The  probability  and  cumulative  density  functions  (PDF  and  CDF)  for  a 
Weibull  distribution  are  given  by,  respectively,  by  the  equations 


(1) 


Isolating  x  in  the  CDF  gives  the  equation 


x  =  — a[Zn(l  —  *)]r 


(2) 
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where  x  is  a  random  number  from  a  Weibull  distribution  and  u  is  a  random  number  from 
a  uniform  distribution.  This  method  is  known  as  the  inverse  transformation  method  and 
is  simple  and  useful  for  generating  random  numbers  from  arbitrary  distributions.  In  this 
notation,  ff  and  y  are  related  to  the  scale  and  shape  parameters  reported  by  Wang  et  al.  in 
the  following  manner: 

vTEl 

a  =  e  *  ,  y  =  fe.  (3) 

Using  the  above  information,  along  with  minimum  and  maximum  nanotube  lengths  of  20 
*****  and  800  *****,  respectively,  the  length  of  each  nanotube  in  the  proposed  model  is 
randomly  generated  from  a  Weibull  distribution  using  the  inverse  transformation  method. 
A  histogram  of  1,000  random  lengths  generated  using  these  Weibull  parameters  is  shown 
in  Figure  1  to  demonstrate  the  effectiveness  of  this  technique. 


Histogram  of  Weibull  Distribution  for  SWCNT  Length 


Nanotube  Length  (nm) 


Figure  1.  A  histogram  of  1,000  Weibull  randomly  generated  nanotube  lengths. 
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Other  authors  have  reported  a  distribution  of  nanotube  diameters  in  addition  to  nanotube 
lengths.  Theoretical  diameters  of  carbon  nanotubes  range  from  about  0.5  Tim  to  7.5  Tim 
as  reported  by  Pipes  et  al.  [55],  Experimentally  fabricated  nanotube  diameters  of  more 
than  5  tsth  have  been  reported,  but  an  average  diameter  between  one  and  two  nanometers 
is  reported  by  several  authors  in  any  given  batch  of  nanotubes,  for  a  variety  of  fabrication 
methods  [56,57],  In  addition  to  length  distribution,  Ziegler  et  al.  [30]  also  reported  a 
distribution  of  individual  nanotube  diameters  which  ranges  from  approximately  0.5  nm. 
to  2  BBS.  The  general  PDF  and  CDF  for  a  lognormal  distribution  are  given,  respectively, 
by  the  equations 

=  **  and  F(x)=  V2  T  V2  erf  (4) 

Solving  the  CDF  for  x  yields  the  equation 

x  —  *rftnvi 2u—  l]  +$[)  ^ 

where  x  is  a  random  number  from  a  lognormal  distribution  and  u  is  a  random  number 
from  a  uniform  distribution.. 

Ziegler  et  al.  did  not  report  the  o  and  u  parameters  for  their  lognormal  calculation,  but 
they  did  report  a  histogram  of  their  diameter  distribution  data.  Using  Matlab’s  “dfittool” 
capability,  a  distribution  fit  is  performed  on  their  histogram  to  calculate  these  parameters. 
This  fit  is  shown  in  Figure  2. 
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Figure  2.  Lognormal  distribution  fit  to  the  diameter  data  reported  by  Ziegler  et  al.  [29]. 

The  parameters  of  u  and  o  are  found  to  be,  respectively,  0.02847  and  0.33637.  Using 
these  parameters,  a  random  diameter  for  each  nanotube  was  generated,  with  a  lower  limit 
of  0.5  nm,  as  supported  by  Ziegler  et  al.  and  by  theoretical  computations.  It  is  important 
to  note,  however,  that  this  diameter  is  not  the  effective  diameter  that  is  incorporated  into 
the  calculation  of  the  volume  that  each  nanotube  occupies.  When  dealing  with 
nanocomposites,  the  actual  distance  of  atomic  bonds  must  be  considered  because  it  is  of 
significant  size  compared  to  the  diameter  of  the  particles.  In  this  particular  case,  the 
carbon-carbon  bond  spacing  between  the  nanotubes  and  the  polymer  matrix  must  be 
considered  in  the  calculation  of  the  effective  volume  occupied  by  the  nanotubes.  A 
carbon-carbon  bond  is  approximately  0.34  mn  in  length,  thereby  increasing  the  effective 
diameter  of  each  nanotube  by  0.68  mn.  A  diagram  of  this  effect  is  shown  in  Figure  3. 
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Figure  3.  The  effective  diameter  of  SWCNT  as  reinforcement  is  increased  due  to  the  carbon-carbon 
bonding  between  the  nanotube  and  the  matrix.  The  actual  nanotube  diameter  is  depicted  as  R_n  and 
the  effective  diameter  including  the  carbon-carbon  spacing  is  shown 


This  effect  is  incorporated  into  the  random  diameter  generation  but  is  ignored  in  the 
random  length  generation.  This  is  done  because  the  0.68  Tan  length  scale  is  significant 
when  compared  to  the  average  diameter  of  approximately  1  mm  but  not  significant  when 
compared  to  the  average  length  of  approximately  200  nm.  A  histogram  of  2,000  random 
effective  diameters  is  shown  in  Figure  4  to  demonstrate  this  distribution. 
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Lognormal  Distribution  of  Effective  Nanotube  Diameters  Including  Carbon-Carbon  Spacing 


2  25  3  35 

Effective  Nanotube  Diameter  in  Matrix  (nm) 


Figure  4.  A  histogram  of  the  lognormal  random  effective  diameters  generated  using  the  parameters 
from  Ziegler’s  paper  [29]  are  presented,  including  carbon-carbon  spacing. 


The  final  parameter  that  is  incorporated  into  the  geometry  generation  of  the  nanotubes  is 
waviness.  As  mentioned  in  the  introduction,  waviness  can  be  considered  a  function  of 
length.  To  incorporate  this  into  the  model,  each  nanotube  was  divided  into  ten  segments, 
which  are  computationally  treated  as  individual  fibers.  Each  segment  is  straight  and  there 
is  an  angle  of  deviation  between  each  segment.  This  angle  is  random,  varying  uniformly 
between - 4^  and  -S^3L  ,  where  ®max  is  some  upper  limit.  The  upper  limit  of  the 

angle  is  determined  by  the  nanotube  length.  Due  to  insufficient  experimental  data  on  the 
subject,  a  linear  model  of  waviness  is  assumed  and  thus  the  function  for  the  maximum 
angle  of  deviation  between  any  two  segments  for  a  particular  nanotube  is  given  as 

0__,  = - — — —  X  actual  nanotube  length.  i  M 

max  uppsrkncth.  limit  &  W 
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With  these  concepts  in  place,  the  geometry  generation  routine  is  implemented  with 
relative  ease.  The  positions  of  each  nanotube  are  generated  randomly  throughout  the 
RVE  and  any  time  a  nanotube  segment  ends  outside  of  the  RVE,  that  nanotube  is  ended 
and  a  new  nanotube  is  begun.  With  the  diameter  and  length  of  each  nanotube  stored,  the 
volume  of  each  nanotube  is  easily  calculated.  The  generation  process  continues  until  the 
desired  volume  fraction  is  reached.  Any  individual  complete  geometry  generation  is 
referred  to  as  a  microstructure.  A  single  RVE  microstructure  with  0.0089%  volume 
fraction  (1%  wt.)  of  nanotubes  in  place  is  illustrated  in  Figure  5. 


Figure  5.  (a)  A  single,  full  RVE  with  0.0889%  v.f.  nanotubes  is  depicted,  (b)  A  close-up  view  of  one 
section  of  the  RVE  is  shown.  The  waviness  of  the  nanotubes  is  easily  observed.  Longer  nanotubes 
are  wavier  than  shorter  nanotubes  to  simulate  the  physical  phenomenon. 

It  should  also  be  noted  that  most  authors  report  the  weight  fraction  of  nanotubes  rather 
than  the  volume  fraction.  The  conversion  between  the  two  quantities  is  not  difficult  as 
long  as  the  densities  of  the  nanotubes  and  matrix  material  are  known.  These  are  input 
parameters  that  are  easily  changed  in  the  model  depending  on  what  type  of  nanotubes  and 
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what  composite  matrix  material  are  desired.  The  conversion  from  weight  fraction  to 
volume  fraction  is  given  by  the  equation  [55] 

tl'm  wfntPmatrix-*'^~wfnt)Pnt 

where  w/nt.  is  the  weight  fraction  of  nanotubes,  Pma?r:x  is  the  density  of  the  polymer 
matrix,  and  is  the  density  of  the  nanotubes. 

2.2  Fiber  Partitioning 

After  generating  the  nanotube  geometry  for  any  single  microstructure,  it  becomes 
necessary  to  partition  each  nanotube  because  of  the  nature  of  the  fiber  embedment. 
Partitioning  is  a  process  in  which  reinforcing  fibers  that  cross  the  boundaries  of  one 
element  into  another  element  are  divided  at  that  boundary.  In  finite  element  analysis,  the 
RVE  is  divided  into  many  elements.  These  elements  allow  for  the  calculation  of 
mechanical  displacements,  temperature  values,  or  any  other  number  of  quantities, 
depending  on  the  fundamental  equation  you  are  seeking  to  solve,  to  be  found  at  finite 
points  throughout  the  RVE.  The  accuracy  of  FEA  increases  with  the  number  of  elements, 
but  each  element  adds  additional  degrees  of  freedom  to  the  system  that  must  be  solved, 
thereby  increasing  the  computational  time  required  to  solve  the  problem.  After  the  initial 
geometry  generation,  each  nanotube  is  divided  into  ten  segments  to  incorporate  waviness, 
as  described  above.  A  visualization  of  this  initial  division  is  shown  in  Figure  6. 
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Figure  6.  The  nanotube  geometry  of  a  small  section  of  an  RVE  before  liber  partitioning  is  presented. 
The  nanotubes  are  divided  into  ten  segments  each,  with  an  “x”  marking  the  end  of  each  segment. 

However,  to  perform  FEA,  only  the  properties  of  the  matrix  and  the  reinforcing  fibers 
that  lie  in  each  particular  element  are  of  use.  Thus  a  fiber  that  begins  in  one  element  and 
ends  in  a  neighboring  element  would  need  to  be  partitioned  at  the  point  it  crosses  the 
element  boundary  into  two,  separate,  shorter  fibers.  If  the  entire  fiber  is  included  in  every 
element  in  which  it  passes,  the  reinforcing  properties  of  that  fiber  would  be  grossly 
overestimated.  It  can  be  imagined  then  that  one  fiber  can  potentially  cross  several 
element  boundaries,  thus  being  divided  several  times,  especially  as  the  number  of 
divisions  within  the  RVE  increases.  Though  this  is  not  physically  precise,  this  method 
has  been  shown  to  be  computationally  accurate  and  must  be  used  for  the  Embedded  Fiber 
Method  (EFM),  which  will  be  discussed  in  the  following  section.  The  partitioning 
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algorithm  used  in  the  current  model  is  taken  from  the  work  of  Esteva  [42]  and  based  upon 
the  work  of  Ranjbaran  [44,45],  It  is  important  to  note  that  although  the  libers  are  divided 
into  essentially  smaller  fibers  by  partitioning,  none  of  the  waviness  and  length 
distribution  information  is  lost  because  this  is  all  taken  to  account  in  the  geometry 
generation  process.  A  sample  microstructure  is  presented  in  Figure  7  to  illustrate  the 
concept  of  fiber  partitioning. 


Figure  7.  (a)  A  visualization  of  the  same  microstructure  as  Fig.  5  with  “x’s”  marking  the  division  of 
individual  nanotubes  into  segments,  (b)  A  magnified  view  of  one  region  of  the  RVE.  The  circles 
represent  segments  that  were  partitioned  because  they  crossed  element  boundaries. 

After  generating  the  geometry  and  partitioning  the  nanotubes  for  an  RVE,  the  final  step 
before  solving  the  system  using  the  Finite  Element  Method  is  to  incorporate  the 
reinforcing  properties  of  the  nanotubes  into  those  of  the  bulk  matrix  material.  This  is 
accomplished  using  the  Embedded  Fiber  Method. 
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2.3  The  Embedded  Fiber  Method 

As  mentioned  in  the  introduction,  the  Embedded  Fiber  Method  (EFM)  was  originally 
developed  for  the  Finite  Element  Analysis  of  steel-rod  reinforced  concrete  by  Ranjbaran 
[44,45],  The  fundamental  principle  behind  this  method  is  the  assumption  of  perfect 
bonding  between  the  bulk  matrix  material  and  the  reinforcing  structures,  thus  allowing 
the  direct  addition  of  the  stiffness  characteristics,  or  any  other  property,  to  those  of  the 
matrix  material.  With  the  recent  advancements  in  nanotube  functionalization  and 
interfacial  bonding  mentioned  in  the  introduction,  this  assumption  is  justified  for  the  case 
of  functionalized  CNT  composites.  To  illustrate  the  general  process,  a  brief  description 
of  the  process  for  the  elastic  case  is  presented. 

The  stiffness  characteristics  of  the  fibers  and  matrix  alone  are  calculated  using  the 
traditional  finite  element  formulation.  The  stiffness  matrix  for  an  isotropic  material 
element  is  given  by  the  equation 


[KVrnot  =  hfA'[BY[D][B]Teda  =  jSnt[By[D][B]T*\jldsdt,  (8) 


where  the  strain  matrix  B,  the  material  matrix  D,  and  the  Jacobian  matrix  J  are  given, 
respectively,  by 


AuUx 

dv  / 


(9) 
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Equations  (9),  (10),  and  (11)  represent  the  matrices  for  the  two-dimensional  plane  stress 
FEA  problem,  which  is  what  the  current  model  is  solving.  In  Equation  (10),  E  is  the 
modulus  of  elasticity  of  the  matrix  material  and  v  is  Poisson’s  ratio.  For  a  four-noded 
quadrilateral  element  with  two  degrees  of  freedom  per  node,  there  are  four  shape 
functions  used  to  map  between  the  local  parametric  coordinates  r  and  s  and  the  global 
coordinates  x  andy.  These  shape  functions  are  seen  in  Equation  (9)  and  are  given  by 


[N}={WtNsMsM4]T, 

(12) 

Nt=±tl-rKls), 

(13) 

#2  =  ;(l+r)(l-s), 

(14) 

W3=j(lir)(Hs), 

(15) 

N4=i(l-r)(l-s), 

(16) 

These  shape  functions  allow  the  mapping  of  local  coordinates  to  global  coordinates  using 
the  equations 


,r(r,s)  =  [jV]  [ri]  a ndyfas)  =  [jV][v]  , 


(17) 
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where  u  and  v  are  the  x  and  y  coordinates  of  the  element  nodes,  respectively.  In  other 
words,  by  using  the  shape  functions  in  Equations  (13)  through  (16) ,  the  x  and  y  positions 
corresponding  to  any  parametric  point  ( r,s )  can  be  found  easily  using  Equation  (17). 
Furthermore,  the  local  derivatives  are  mapped  to  global  derivatives  using  the  Jacobian 
from  Equation  (11)  and  the  derivates  of  the  shape  functions.  This  mapping  is  given  by 
the  equation 


(18) 


With  these  derivatives  and  coordinates  in  hand,  the  strain  matrix  and  stiffness  matrix  for 
an  element  are  easily  calculated.  Due  to  the  discrete  nature  of  FEA,  the  integral  of 
Equation  (8)  is  determined  numerically  using  Gaussian  quadrature.  For  any  particular 
element,  the  shape  functions  and  the  Jacobian  are  evaluated  at  each  of  the  Gaussian 
quadrature  points  within  the  element  and  summed.  For  a  two-dimensional  square 
element,  as  is  the  case  in  the  present  model,  these  points  are  given  as 


(19) 


The  stiffness  matrix  can  now  be  written  as 


(20) 


These  finite  element  formulations  can  be  found  in  any  standard  Finite  Element  Method 
book,  such  as  those  by  Akin  [58],  Desai  and  Abel  [59],  and  Zienkiewicz  et  al.  [60]. 

With  the  ability  to  calculate  the  stiffness  matrix  for  the  bulk  matrix  material,  it  is  now 
necessary  to  calculate  the  stiffness  matrices  for  the  reinforcing  fibers.  In  the  current 
model,  the  nanotubes  are  approximated  as  line  elements.  This  approximation  is  accurate 
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and  economical  as  shown  by  Konrad  and  Graovac  [43],  With  this  approximation,  the 
stiffness  matrix  of  a  one-dimensional  fiber,  expanded  using  zeros  to  two  dimensions,  is 
given  as 


ft'  -  — 
hfV>L  ~  L 


1  0-1  0 
0  0  0  0 

-10  10 
.  0  0  0  0. 


(21) 


This  general  stiffness  calculation  is  for  a  fiber  that  is  oriented  in  the  local,  positive  x- 
direction.  In  Equation  (21),  A  is  the  cross-sectional  area  of  the  fiber  and  L  is  the  fiber 
length.  To  transform  this  matrix  from  its  local  system  to  the  global  coordinate  system, 
accounting  for  the  true  alignment  of  the  nanotube,  a  matrix  of  direction  cosines  is  easily 
employed  using  the  equation 

•r«f  n  ft 

(22) 


[CJ= 


cost  sint  0  0 

—sint  cost  0  0 

0  0  cost  sint 

0  0  —sint  cost . 


In  Equation  (22),  t  is  the  angle  between  the  fiber  and  the  positive  x-axis  of  the  global 
coordinate  system. 

The  next  step  to  prepare  the  fiber’s  stiffness  matrix  to  be  embedded  with  the  bulk 
material’s  stiffness  matrix  for  a  single  element  is  to  map  the  physical  coordinates  of  the 
fiber  endpoints  to  the  corresponding  parametric  coordinates.  The  parametric  coordinates 
of  the  fiber  endpoints  are  found  using  a  Newton-Rhapson  scheme  which  is  shown  in 
Appendix  A. 

Finally,  the  full,  global  stiffness  matrix  for  a  single  fiber  in  a  single  element  can  be  found 
by  interpolating  the  fiber’s  properties  that  are  given  in  Equation  (22)  to  the  nodes  of  the 
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element  under  consideration.  This  is  accomplished  by  utilizing  the  transformation  matrix 
T  with  the  equation 


Kfiil  =  [T][A>i2,J[r]T  ,  where 

^  ^2^* i*^i)  ®  ^  ^ 

o  %&!*%>  0  0  fly***i>  0  %(*i**i) 

0  ^(r^Sj)  0  AT3(rSJss)  0  A4(rs,  Sj)  0 

0  ATt(r^s2)  0  iV2(r2j  sv)  d  iV3(r^s:2)  6  %(?2>£2). 


(23) 


In  Equation  (23),  r  and  5  represent  the  local  parametric  coordinates  of  the  fiber  endpoints. 
Having  introduced  this  final  transformation,  the  fiber’s  properties  can  now  be 
incorporated  into  those  of  the  bulk  matrix  material  using  simple  addition: 

[M  [A'] 

771  Gw  +  M/».  (24) 

It  is  important  to  note  here  that  while  the  above  steps  are  described  for  the  case  an  elastic 
problem,  this  method  is  truly  versatile  and  may  be  used  with  minimal  change  for  thermal, 
electrical,  or  various  other  problems. 

The  Embedded  Fiber  Method  is  a  bit  daunting  at  first  glance.  Not  only  are  there 
numerous  equations  to  take  into  consideration,  but  this  process  must  be  accomplished  for 
every  fiber  in  an  element  as  well  as  for  every  element  in  the  RVE.  However,  the  benefits 
achieved  by  accomplishing  this  extra  work  are  indeed  significant.  The  primary 
advantage  is  that  the  EFM  allows  the  use  of  constant  quadrilateral  elements.  Most 
models  that  use  FEA  must  utilize  intense  mesh  refinement  along  the  nanotube/matrix 
boundaries  in  order  to  accurately  capture  the  load  transfer  properties  of  the  interface.  For 
example,  a  picture  of  the  finite  element  mesh  used  by  Fi  and  Chou  [19]  in  their  model  is 
shown  in  Figure  8. 
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Figure  8.  The  finite  element  mesh  of  the  carbon  nanotube/fiber  cross-ply  composite  used  by 
Li  and  Chou  [19].  The  intense  mesh  refinement  along  the  nanotube  matrix  boundaries  is  to 
be  noted. 


By  utilizing  a  constant  quadrilateral  mesh,  accurate  results  can  be  found  without  the 
costly  human  intervention  of  most  other  models.  For  example,  the  model  of  Li  and  Chou 
shown  above  uses  approximately  357,000  elements  and  360,000  nodes  to  achieve  their 
results  [19].  The  fiber-matrix  interfaces  must  be  identified  and  the  mesh  must  be  refined 
in  those  areas.  When  compared  to  the  3,600  elements  and  3,721  nodes  needed  for 
accurate  results  in  the  current  model,  the  advantages  of  the  Embedded  Fiber  Method 
become  quite  clear.  For  the  thermal  model  with  one  degree  of  freedom  per  node,  there 
are  approximately  lOOx  fewer  equations  to  solve  in  the  proposed  model.  There  are  about 
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200x  fewer  equations  to  solve  for  the  elastic  problem.  No  attempt  to  quantify  the  time 
saved  by  eliminating  the  human  intervention  process  is  made  herein.  However,  the 
complete  automation  of  the  process  by  using  the  EFM  is  undoubtedly  useful. 

Without  the  Embedded  Fiber  Method,  it  is  extremely  difficult  to  accurately  model  any 
significant  volume  fraction  of  nanotubes  because  the  computational  time  is  simply  too 
great.  When  considering  non-linear  properties,  the  number  of  equations  to  be  solved  is 
multiplied  by  the  number  of  iterations  or  incremental  steps  employed.  The  Embedded 
Fiber  Method  is  the  only  practical  solution  to  evaluate  the  non-linear  properties  of  a 
CNT-reinforced  composite  on  a  scale  larger  than  only  a  few  nanotubes  and  is  thus 
utilized  by  the  model  proposed  in  this  paper. 
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Chapter  3 


The  Non-Linear  Finite  Element  Model 


Having  established  the  Embedded  Fiber  Method  as  the  procedure  for  incorporating  the 
physical  properties  of  the  nanotubes  into  the  bulk  composite  material  properties,  the 
proposed  method  can  now  be  detailed.  Esteva  [42]  developed  this  model  using  the 
Embedded  Fiber  Method  with  the  goal  to  determine  the  effective  elastic  and  thermal 
properties  of  the  composite.  Using  a  linear  model,  the  RVE  was  displaced  to  5%  strain 
and  the  reactions  at  the  fixed  bottom  nodes  were  used  to  calculate  the  effective  modulus 
of  elasticity.  Similarly,  a  temperature  gradient  was  applied  between  the  upper  and  lower 
boundaries  and  heat  flux  at  the  lower  nodes  was  used  to  calculate  the  effective  thermal 
conductivity.  While  this  model  made  a  significant  advancement  in  the  area  of  composite 
modeling,  one  of  the  primary  tools  used  by  researchers  to  describe  a  new  or  existing 
material  is  its  stress-strain  curve.  By  adopting  a  linear  solution,  it  is  not  possible  to 
accurately  capture  the  non-linear  properties  of  the  polymer  matrix  or  those  of  the 
nanotubes  themselves.  The  proposed  model  explores  both  of  these  non-linearities  and 
how  they  affect  the  thermal  and  elastic  properties  of  the  composite  as  a  whole.  First,  the 
relevant  non-linearities  are  identified.  After  identifying  which  non-linearities  to  consider, 
the  various  non-linear  solution  techniques  are  discussed  and  how  they  are  applied  to  the 
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proposed  model.  Finally,  the  specific  details  of  the  model  are  discussed  along  with  the 
parameters  chosen  for  experimental  comparison. 

3.1  Non-Linearity  of  Carbon  Nanotubes  and  Epoxies 


The  non-linear  elastic  and  thermal  properties  of  carbon  nanotubes  are  well  documented  in 
the  literature,  both  theoretically  and  experimentally.  For  the  non-linear  elastic  properties, 
Tserpes  et  al.  [34]  presented  a  detailed  finite  element  model  of  lone  nanotubes  using  a 
structural  mechanics  approach.  They  identify  the  stress-strain  relationship  for  armchair 
and  zigzag  SWCNTs.  Additionally,  Tserpes  et  al.  [6]  presented  a  progressive  fracture 
model  to  calculate  the  stress-strain  relationship  for  SWCNT’s  with  perfect  structures  and 
for  those  with  imperfections.  Belytschko  et  al.  [7]  developed  a  molecular  mechanics 
model  to  simulate  nanotube  fracture.  They  reported  the  stress-strain  curve  for  perfect  and 
imperfect  zigzag  nanotubes  as  well  as  a  curve  for  various  different  perfect  nanotube 
structures.  Natsuki  et  al.  [4]  described  a  structural  mechanics  approach  to  modeling 
nanotube  elastic  properties.  Finally,  Meo  and  Rossi  [5]  presented  a  molecular-mechanics 
based  FEM  for  the  prediction  of  SWCNT  elastic  properties.  Each  of  these  models  takes  a 
slightly  different  approach,  but  they  all  identify  ~  10-15%  strain  as  the  region  where 
significant  non-linear  behavior  is  first  exhibited  by  SWCNTs. 

The  non-linear  thermal  properties  of  SWCNTs  have  also  been  modeled,  though  less 
frequently  than  the  elastic  properties.  Grujicic  et  al.  [15]  presented  a  molecular  dynamics 
based  model  of  the  thermal  conductivity  of  SWCNTs  in  which  they  investigate  the  non¬ 
linear  dependence  of  thermal  conductivity  on  temperature.  The  strong  dependence  they 
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report  is  depicted  in  Figure  9  and  is  the  motivation  for  the  non-linear  thermal  analysis 
performed  in  this  paper. 


Figure  9.  The  non-linear  thermal  properties  of  SWCNTs  are  shown  in  the  analysis  of  Grujicic  et  al. 
[15]. 

Kawamura  et  al.  [61]  also  reported  the  temperature  dependence  of  SWCNT  thermal 
conductivity  using  a  molecular  dynamics  simulation.  Additionally,  Xu  et  al.  [47] 
reported  this  effect  indirectly  by  reporting  the  temperature  dependence  of  the  thermal 
conductivity  of  composites  reinforced  with  SWCNTs. 

Along  with  the  theoretical  modeling  evidence  presented  above,  experimental  evidence 
corroborates  the  non-linear  thermal  and  elastic  properties  of  SWCNTS.  Yu  et  al.  [8,27] 
performed  tensile  load  testing  using  an  AFM  tip  to  document  the  non-linear  stress-strain 
relationship  of  SWCNTs.  Their  results  show  the  onset  of  non-linear  behavior  to  occur  at 
greater  than  10%  strain.  Walters  et  al.  [26]  also  reported  the  lower  bound  of  the  yield 
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strain  of  SWCNT  ropes  to  be  about  6%.  For  thermal  properties,  Shaikh  et  al.  [62] 
reported  a  range  of  thermal  conductivities  for  CNT  ropes  to  range  from  0.1-2,000  — — . 

Epoxy  composites  exhibit  non-linear  elastic  properties  at  a  much  lower  strain  level.  Zhu 
et  al.  [17,22]  present  stress-strain  curves  for  neat  epoxy  in  which  non-linearity  is  shown 
at  3%  strain.  Li  et  al.  [63]  present  similar  effects  with  a  load-displacement  curve  of  neat 
epoxy.  Sun  et  al.  [46]  and  Gojny  et  al.  [54]  also  indicate  3%  strain  as  the  onset  of  non¬ 
linearity  in  neat  epoxy.  Each  author  indicates  epoxy  failure  at  less  than  6%  without 
reinforcement.  The  significantly  lower  onset  of  non-linear  behavior  and  failure  in  the 
epoxy  than  in  the  lone  nanotubes  led  to  the  conclusion  that  the  non-linear  behavior  of  the 
nanotubes  may  be  ignored  and  that  of  the  epoxy  must  be  taken  into  account  for  an 
accurate  elastic  model. 

In  general,  polymers  are  reported  to  have  constant  values  of  thermal  conductivity.  Xu  et 
al.  [47]  reported  a  constant  thermal  conductivity  of  ~  0.2  _..for  their  PVDF  polymer. 

Wang  et  al.  [48]  reported  a  constant  thermal  conductivity  of  ~  0.18  — -  for  the  neat 
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epoxy  resin  used  in  their  experiments.  Moisala  et  al.  [16]  reported  a  constant 
~  0.255  — —  for  epoxy  consisting  of  a  bisphenol-A  resin  and  an  aromatic  amine  hardener. 

Du  et  al.  [14]  also  reported  a  constant  thermal  conductivity  of  ~.18  — -  for  their  PMMA 
polymer.  Finally,  Cai  and  Song  [49]  used  a  semi-crystalline  PU  dispersion  with  a 
constant  thermal  conductivity  of  ~  0.15  — — .  Because  of  the  lack  of  any  indication  of  a 

non-linear  temperature  dependence  of  the  thermal  conductivity  of  the  neat  epoxy,  only 
the  non-linear  thermal  conductivity  of  the  nanotubes  should  be  considered  for  a  thermal 
model  of  CNT-reinforced  polymer  composites. 
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Based  upon  this  review  of  the  non-linear  thermal  and  elastic  properties  of  both  carbon 
nanotubes  and  polymer  matrix  materials,  it  is  concluded  that  the  non-linear  elastic 
properties  of  the  epoxy  matrix  will  be  included  in  the  elastic  analysis  of  the  model 
presented  in  this  paper  and  the  non-linear  thermal  properties  of  the  CNTS  will  be 
included  for  the  thermal  analysis.  With  these  choices  made,  the  next  step  is  to  select  the 
non-linear  approach  to  take  in  each  case. 


3.2  Non-Linear  FEA  Solution  Techniques 


In  general,  there  are  two  main  sources  of  non-linearity  in  finite  element  analysis,  as 
described  in  references  such  as  Desai  and  Abel  [59].  The  first  type  is  called  material  or 
physical  non-linearity.  This  type  of  non-linearity  stems  from  constitutive  laws  and  is 
inherent  in  the  physical  properties  of  the  material.  The  second  kind  of  non-linearity 
relates  to  geometrical  considerations.  When  materials  undergo  large  deformations,  the 
geometry  of  the  material  may  change  drastically  and  this  can  alter  the  equations  used  to 
govern  the  finite  element  analysis  of  that  body  [59].  In  both  the  thermal  and  elastic  cases 
considered  in  the  present  model,  material  non-linearity  is  the  only  aspect  considered. 

For  material  non-linearity,  Desai  and  Abel  [59]  outline  three  kinds  of  non-linear  solution 
techniques  commonly  used  for  finite  element  analysis:  the  incremental  approach,  the 
iterative  approach,  and  a  mixed  approach.  In  the  incremental,  or  step-wise,  approach,  the 
load  applied  to  the  system  is  subdivided  into  many  small  increments.  In  general,  these 
increments  do  not  need  to  be  of  the  same  size.  With  the  application  of  each  load 
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increment,  all  the  constitutive  laws  pertaining  to  the  system  are  assumed  to  be  linear. 
The  results  obtained  after  each  step  are  then  used  to  update  the  values  of  the  constants 
used  in  the  previous  linear  increment  based  upon  a  known  relationship.  For  a  non- 
homogeneous  material,  the  constants  may  vary  throughout  the  body.  When  applying  the 
iterative  approach,  the  structure  is  fully  loaded  during  each  iteration.  After  solving  the 
fully  loaded  system,  the  constants  are  then  updated  based  upon  the  results.  In  general, 
equilibrium  is  not  satisfied  after  each  iteration.  The  non-balanced  load  is  then  used  to 
correct  the  solution  until  equilibrium  is  achieved  or  is  sufficiently  close  based  upon 
predetermined  criteria.  Finally,  a  mixed  approach  is  combination  of  the  iterative  and 
incremental  approaches.  Generally,  incremental  loads  are  applied  and  a  number  of 
iterations  are  applied  to  each  load  step  [59]. 

After  careful  consideration  of  the  above  methods,  the  incremental  technique  was  chosen 
for  the  elastic  analysis  to  be  performed  and  both  techniques  were  studied  for  the  thermal 
analysis.  The  incremental  technique  was  chosen  both  for  its  versatility  and  because  it 
most  closely  resembles  the  physical  phenomenon  of  a  tensile  strain  test.  The  iterative 
technique  is  not  applicable  in  the  elastic  analysis  of  the  current  model  because  it  does  not 
provide  any  information  about  the  stresses  and  strains  at  intermediate  load  increments. 
Thus,  a  stress-strain  relationship  cannot  be  obtained  from  an  iterative  approach  in  the 
elastic  case.  In  addition,  the  availability  of  experimental  data  for  comparison  contributed 
to  this  choice. 

For  the  thermal  analysis,  both  an  incremental  approach  and  an  iterative  approach  were 
taken  for  comparison.  An  iterative  approach  is  more  readily  implemented  than  an 
incremental  approach  and  converges  more  rapidly  than  the  incremental  approach.  For 
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this  reason,  it  was  utilized  to  calculate  an  effective  thermal  conductivity  for  the 
composite.  However,  there  is  no  guarantee  that  an  iterative  technique  will  converge  and 
it  provides  no  conductivity  information  at  intermediate  temperature  increments,  so  an 
incremental  approach  was  also  used. 

Upon  deciding  upon  the  approaches  to  select  for  the  thermal  and  elasticity  analyses,  a 
non-linear  finite  element  model  was  developed,  using  the  Embedded  Fiber  Method  as  a 
foundation.  After  generating  the  nanotube  geometry,  partitioning  the  fibers,  and 
embedding  the  fibers,  the  stiffness  matrix  of  the  entire  RVE  is  assembled  using  the 
standard  finite  element  procedures  outlined  in  Chapter  2.  The  standard  force  balance 
equation  for  FEA  is  given  by  the  equation 

[A'Jw  =  F,  (25) 

where  [JT]  is  the  system  stiffness  matrix  assembled  from  the  element  stiffness  matrices 

•  •  •—  ft  #  # 

defined  in  Equation  (24),  m  is  the  displacement  or  temperature  vector  for  all  of  the  mesh 
nodes  in  the  RVE,  and  F  is  the  reaction  force  or  temperature  gradient  vector  at  the  nodes. 
For  the  elasticity  analysis  case,  the  general  finite  element  problem  can  be  solved  in  one  of 
two  standard  ways.  The  first  approach  is  to  prescribe  a  force  to  be  applied  to  the 
structure  and  then  to  calculate  the  resulting  displacements  at  all  the  nodes.  The  second 
approach  is  to  prescribe  displacements  to  specific  nodes  and  then  calculate  the  resulting 
displacements  of  all  the  other  nodes.  In  each  method,  the  final  displacements  can  be  used 
to  recover  the  reaction  forces  at  each  of  the  fixed  nodes.  The  current  model  uses  the 
displacement-controlled  approach  because  the  uniform  displacement  of  an  entire  edge  is 
needed  to  generate  a  stress-strain  curve  for  the  RVE  as  a  whole. 
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To  implement  the  incremental  solution  technique,  displacement  increments  were  applied 
to  the  elasticity,  problem  and  temperature  increments  were  applied  to  the  thermal 
problem.  The  number  of  increments  necessary  for  an  accurate  solution  is  one  of  the 
results  of  the  present  report.  After  each  increment  is  applied,  the  resulting  displacements 
or  temperatures  of  every  node  in  the  mesh  are  calculated. 

For  the  elasticity  case,  there  are  two  degrees  of  freedom  per  node;  x  and  y  displacements. 
With  the  displacements  at  each  node  in  hand  after  an  incremental  displacement  is  applied, 
the  strain  at  each  node  can  be  calculated  and  interpolated  to  the  Gaussian  quadrature 
points  using  the  strain  matrix  B  from  Equation  (9),  and  the  equation 

*  =  [£]«.  (26) 
In  two-dimensions,  the  strain  vector  *  is  a  three  component  vector  including  the  strain  in 
the  x  and  y  directions,  sx  and  as  well  as  the  shear  stain,  .  With  the  strain,  the  stress 
vector  at  each  Gauss  point  can  be  calculated  using  the  relation 

=  Ml*.  (27) 

Similar  to  the  strain  vector,  the  stress  vector  a  in  two-dimensions  is  a  three  component 
vector  including  the  strain  in  the  x  and  y  directions,  &x.  and  as  well  as  the  shear  stress, 
T™.  Recalling  the  earlier  description  of  the  basic  finite  element  method,  [13]  is  the 
material  matrix  given  in  Equation  (10).  It  is  worth  pausing  here  to  describe  the  modulus 
of  elasticity  E  that  appears  in  the  material  matrix  and  how  it  is  calculated  in  the  current 
model.  In  any  linear  model,  the  modulus  of  elasticity  of  a  solid  material  is  assumed  to  be 
constant.  For  a  perfectly  elastic  material,  this  modulus  can  be  calculated  as  the  slope  of 
the  stress-strain  curve  for  that  material.  Though  no  material  is  perfectly  elastic,  this 
assumption  is  valid  for  most  “nearly”  elastic  materials.  For  materials  that  exhibit  non- 
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linearity,  the  slope  of  the  stress-strain  curve  changes  as  stress  levels  change.  For 
materials  that  exhibit  plasticity,  the  curve  decreases  at  the  onset  of  plastic  deformation. 
For  work-hardening  materials,  the  slope  of  the  curve  may  even  increase  at  higher  stress 
or  strain  levels.  The  slope  of  the  curve  at  any  given  strain  level  is  known  as  the  tangent 
modulus.  In  this  report,  the  tangent  modulus  of  a  polymer  material  is  a  function  of  strain, 
not  assumed  to  be  constant,  and  will  simply  be  referred  to  as  the  modulus  or  elastic 
modulus  of  the  material. 

In  general  plasticity  theory,  there  are  two  competing  ways  to  categorize  whether  a 
material  has  begun  to  undergo  plastic  deformation.  The  Tresca  and  von  Mises  criteria  are 
two  different  methods  which  take  into  account  the  full  stress  state  at  one  specific  point  of 
a  material  to  come  up  with  a  single  value  which  describes  the  effective  stress  at  that 
point.  The  von  Mises  criterion  is  the  more  conservative  of  the  two  and  is  thus  employed 
in  the  current  model.  The  von  Mises  equivalent  stress  in  two  dimensions  is  given  by  the 
equation  [59,64] 

°V'Arr  =  JV 2  [fe  “  ayT  +  ay  +  ^  +  6Tr/]  (28) 

After  finding  the  strain  vector  at  a  node  using  Equation  (26)  and  using  it  to  find  its 
corresponding  stress  vector  by  Equation  (27),  the  von  Mises  stress  at  each  Gauss  point  in 
each  element  can  be  calculated  using  Equation  (28).  The  von  Mises  stress  is  then  used  to 
refer  to  the  stress-strain  curve  of  the  neat  epoxy  to  calculate  the  tangent  modulus  at  any 
particular  Gauss  point.  In  this  manner,  the  modulus  at  each  Gauss  point  throughout  the 
RVE  can  be  updated  after  each  displacement  increment.  This  procedure  allows  the  non¬ 
linearity  of  the  material  to  be  effectively  incorporated  into  the  model.  A  synopsis  of  the 
incremental  procedure  is  presented  in  Figure  10. 
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<7 

Update  Elastic  Modulus  Using  Material  Stress-Strain  Curve  at  Each  Gauss  Point 

Calculate  Stiffness  Matrixfor  Element 

Update  Stresses  at  each  Gauss  Point  using  Interpolated  Displacements 

<7 

Add  Nanotube  Stiffness  Matrix  to  Material  Stiffness  Matrixfor  That  Element 

Assemble  System  Stiffness  Matrix 

Apply  Boundary  Conditions 

Solve  System  for  Nodal  Displacements 

U 

End:  Post  Processing 

Figure  10.  A  block  diagram  of  the  non-linear  incremental  scheme  for  the  elastic  case  is  presented. 
Notice  the  two  major  loops  which  are  the  driving  force  behind  the  computational  time  required  by 
the  model. 

For  the  thermal  problem,  the  incremental  approach  follows  essentially  the  same  approach 
as  the  elastic  problem,  with  only  a  few  differences.  The  first  difference  between  the 
thermal  and  elastic  problem  is  that  the  thermal  case  requires  only  one  degree  of  freedom 
per  node,  temperature.  This  in  turn  greatly  reduces  the  computational  cost  to  solve  the 
problem  because  the  number  of  equations  to  solve  is  halved.  The  second  difference  is 
that  instead  of  updating  the  elastic  modulus  of  the  material,  the  thermal  conductivity  of 
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the  nanotubes  is  updated  after  each  incremental  temperature  increase.  This  process  is 
simply  performed  with  the  knowledge  of  the  temperature  at  each  of  the  nodes  throughout 
the  RVE.  Using  the  shape  functions  given  in  Equations  (13)  through  (16),  the 
temperature  is  interpolated  from  the  nodal  temperatures  to  define  a  temperature  at  each  of 
the  nanotube  endpoints.  Because  the  nanotube  thermal  conductivity  is  quite  high, 
approximately  four  orders  of  magnitude  greater  than  that  of  the  polymer  matrix,  an 
average  of  the  two  nanotube  endpoint  temperature  is  assumed  as  the  constant  temperature 
along  that  nanotube.  The  thermal  conductivity  of  that  particular  nanotube  is  then  updated 
based  upon  the  curve  of  Grujicic  et  al.  [15]  shown  in  Figure  9.  The  thermal  conductivity 
is  updated  using  the  equation 

Knt  =  3E~ST5  +  7E~ST 4  -  0.0011  T3  +  0.1568712  -  26.5037  +  2471.7,  (29) 

where  Ent  is  the  thermal  conductivity  of  the  nanotube  and  T  is  the  temperature  of  the 
nanotube  in  degrees  Celsius.  This  process  is  repeated  for  every  nanotube  in  each  element 
and  is  repeated  after  each  incremental  temperature  step. 

In  addition  to  the  incremental  approach,  the  iterative  technique  is  also  applied  to  the 
thermal  problem.  In  this  process,  the  finite  element  solution  is  solved  in  one  total 
temperature  increase.  Using  the  temperature  data  recovered  from  this  solution,  the 
thermal  conductivities  are  updated  using  the  same  procedure  as  the  incremental  approach. 
The  problem  is  then  solved  again  using  this  data.  This  process  is  repeated  until  the 
solution  converges  to  a  solution.  The  number  of  iterations  required  for  convergence  is 
one  of  the  results  presented  in  the  discussion  section  of  the  paper. 

Overall,  there  are  two  aims  to  the  non-linear  analysis  that  is  performed  by  the  proposed 
model.  The  first  aim  is  to  create  a  stress-strain  curve  for  any  CNT -reinforced  polymer 
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composite.  The  second  aim  is  to  discover  whether  or  not  the  non-linear  thermal 
properties  of  carbon  nanotubes  need  to  be  included  in  thermal  models  of  CNT/polymer 
composites.  The  results  of  incremental  and  iterative  non-linear  approaches  are  compared 
for  this  aim  to  see  if  one  is  superior  to  the  other.  However,  whatever  the  results  of  the 
model  may  be,  the  success  of  these  aims  is  ultimately  determined  by  whether  or  not  the 
model  accurately  represents  available  experimental  results.  Due  to  this  fact,  the  present 
model  is  compared  to  several  experiments  by  incorporating  the  parameters  reported  in 
those  studies  into  the  details  of  the  model. 

3.3  Details  of  the  Model 

With  the  general  procedure  detailed  above,  the  two  significant  attributes  of  the  proposed 
model  should  be  evident.  First  and  foremost,  the  model  offers  the  option  to  to 
incorporate  the  non-linear  properties  of  both  the  carbon  nanotubes  and  the  polymer 
matrix  into  both  thermal  and  elasticity  finite  element  analysis.  Though  only  the 
polymer’s  elastic  non-linearity  and  the  nanotubes’  thermal  non-linearity  are  incorporated 
into  the  present  elasticity  and  thermal  analyses,  respectively,  it  would  take  only  minimal 
modification  to  incorporate  both  the  non-linear  elastic  properties  of  the  nanotubes  and  the 
thermal  properties  of  the  matrix  should  the  specific  physical  situation  require  it.  For 
instance,  if  a  different  polymer  or  other  material  were  used  in  which  non-linearity  was 
initiated  at  a  higher  strain  level,  this  would  make  the  non-linear  elastic  properties  of  the 
nanotubes  relevant,  and  the  same  procedure  utilized  for  the  present  thermal  analysis 
could  quickly  be  used  for  the  new  elasticity  problem.  The  same  applies  for  a  new 
thermal  situation.  The  second  significant  attribute  of  the  proposed  model  is  its  inherent 
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versatility.  If  it  is  desired  to  model  a  composite  fabricated  with  short-straight  nanotubes, 
such  as  the  ones  developed  by  Li  et  al.  [65]  or  Jia  et  al.  [66],  all  that  is  needed  is  a  change 
in  the  parameters  of  the  length  distribution  described  in  Chapter  2,  Equation  (1).  This 
will  in  turn  make  the  waviness  insignificant  and  short,  straight  nanotubes  will  be 
generated.  If  it  is  desired  to  model  a  composite  with  aligned  nanotubes,  all  that  is  needed 
is  a  minor  change  in  the  geometry  generation  section  of  the  model  in  which  the  first 
nano  tube  segment  is  always  generated  in  the  same  direction.  Finally,  the  model  is 
flexible  to  almost  any  physical  difference  in  properties  among  the  materials  used.  The 
stress-strain  curve  used  to  update  the  elastic  modulus  of  the  matrix  material  is  a  user- 
defined  parameter  and  can  be  readily  changed  depending  upon  the  material  used.  If  a 
new  or  different  curve  for  the  thermal  conductivity  variation  of  SWCNTs  with 
temperature  is  specified,  the  equation  used  to  update  the  thermal  conductivity  of  the 
nanotubes  readily  altered. 

With  these  capabilities  in  mind,  the  elastic  results  of  the  model  are  compared  to  two 
different  experimentally  determined  stress-strain  curves.  The  first  experiment  is  the 
tensile  test  performed  by  Sun  et  al.  [46]  on  an  epoxy  composite  reinforced  in  several 
different  ways.  The  epoxy  used  was  the  EPIK-OTE™  862  resin  cured  with  the 
EPIKURE™  W  curing  agent.  The  stress-strain  curves  for  the  neat  epoxy,  epoxy  with 
pristine  nanotubes,  and  epoxy  with  functionalized  nanotubes  are  reported  and  shown  in 
Figure  11. 
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Figure  11.  Tensile  stress  versus  strain  curves  for  the  various  composites  reported  by  Sun  et  al.  are 
presented  |45|.  The  bottom,  black  line  is  used  in  the  current  model  to  update  the  elastic  modulus  of 
the  epoxy  matrix. 


The  bottom,  black  curve  for  the  neat  epoxy  shown  in  Figure  1 1  is  used  as  the  reference 
curve  by  which  the  modulus  of  elasticity  is  updated  in  the  incremental  elastic  analysis  of 
the  current  model.  One  of  the  fundamental  assumptions  of  the  Embedded  Fiber  Method 
is  that  the  reinforcing  fibers  are  perfectly  bonded  to  composite  matrix,  as  explained  in 
Chapter  Two.  Due  to  this  assumption,  the  results  of  the  model  will  be  compared  to  the 
top,  blue  stress-strain  curve  for  the  epoxy  reinforced  with  1  wt%  functionalized 
nano  tubes. 

Using  a  coordinate  finding  program  developed  at  Rice  University  [42],  the  data  points  of 
the  neat  epoxy  curve  in  Figure  1 1  were  extracted  from  the  reported  figure.  Using  those 
data  points  and  Matlab’s  built-in  polynomial  fitting  routine,  “poly fit”,  the  cubic  stress- 
strain  polynomial  is  given  by  the  equation 
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a  =  -ZBHE^e^  -  1S73M10e2  +  2.864 E*s,  (30) 

with  an  fit  value  of  0.999.  Taking  the  derivative  of  this  function,  the  instantaneous 
tangent  modulus  of  the  epoxy  matrix  is  given  by  the  equation 

^7  =  Eta*  =  —1.193  E12s2  —  3.148E10s  +  2.864E9.  (31) 

The  elasticity  problem  results  of  the  proposed  model  are  also  compared  to  the  work  of 
Zhu  et  al.  [22],  In  their  experiment,  they  used  a  similar  epoxy  and  curing  agent  to  Sun  et 
al.  [46]  and  they  also  used  1  wt%  SWCNT  loading.  However,  they  report  a  different 
stress-strain  relationship  for  the  neat  epoxy,  allowing  the  simulation  to  be  continued  out 
to  5%  strain.  This  curve  is  shown  in  Figure  12. 


Figure  12.  Tensile  stress  versus  strain  curves  for  the  composites  reported  by  Zhu  et  al.  are  presented 
[22].  The  middle,  neat  epoxy  line  is  used  in  the  current  model  to  update  the  elastic  modulus  of  the 
epoxy  matrix. 
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Using  the  epoxy  curve  that  Zhu  et  al.  report,  the  coordinate  finding  program  of  Esteva 
[42]  was  again  employed  and  the  data  points  for  that  curve  were  extracted.  Again  using 
Matlab’s  “polyfit”  routine,  the  stress-strain  polynomial  fit  to  the  data  was  found  to  be 
given  by  the  equation 

a  =  — 1.1084E11s3  -  4.0Ess2  -f  1.93 E9e.  (32) 

This  stress-strain  relationship  yielded  a  tangent  modulus  equation  given  by 

—  =  Eccn  =  -332SE11e2  -  8.0 E8e  +  1.93F9.  (33) 

It  should  be  noted  here  that  in  order  to  calculate  the  tangent  moduli  in  Equations  (31)  and 
(33),  a  single  value  of  strain  is  needed.  The  calculation  of  a  single  stress  was  explained 
in  the  above  section  utilizing  the  von  Mises  criterion,  but  no  mention  was  made  of 
calculating  a  single  strain.  There  are  a  few  references  to  a  von  Mises  strain,  but  this 
quantity  is  almost  never  used  and  several  different  reported  equations  for  the  von  Mises 
strain  were  found.  Taking  this  into  account,  the  most  direct  solution  to  the  problem  was 
to  calculate  the  von  Mises  stress  using  Equation  (28)  and  then  to  solve  Equation  (30)  for 
the  strain.  In  this  manner,  one  could  use  the  von  Mises  stress  to  calculate  a  single  strain 
value.  To  solve  Equation  (30)  for  the  strain,  a  simple  Newton-Rhapson  root-finding 
algorithm  was  developed  to  find  the  strain  for  any  given  von  Mises  stress.  This  value  is 
then  used  in  Equation  (31)  to  find  the  tangent  modulus. 

It  is  also  worth  mentioning  here  that  an  additional  limitation  was  placed  upon  this  routine 
to  ensure  success.  As  noted  in  the  composite  materials  literature,  it  is  possible  and 
probable  that  local  stress  and  strain  concentrations  will  occur  throughout  the  RVE  due  to 
the  local  reinforcement  provided  by  the  nanotubes  [19,67].  In  reality,  when  the  strain 
becomes  too  high,  the  reinforcing  agent  will  tear  from  the  reinforcing  matrix  and 


45 


reinforcement  will  drastically  decrease.  This  effect  can  be  seen  visually  in  the 
experimental  work  of  Sun  et  al.  [46],  Zhu  et  al.  [22],  and  Gojny  et  al.  [54],  but  is  most 
readily  observed  in  the  work  of  Li  and  Chou  [19]  and  shown  in  Figure  13. 


Figure  13.  A  close  view  of  the  stress-concentrations  at  nanotube  tips  due  tensile  strain  is  presented. 

Li  and  Chou  developed  this  picture  from  their  damage  sensing  model  |19). 

Since  the  proposed  model  relies  upon  the  stress-strain  curve  reported  by  experimentalists 
to  determine  the  modulus  of  elasticity,  values  of  strain  beyond  the  upper  limit  of  the 
reported  stress  strain  curve  present  a  problem.  The  quadratic  polynomial  of  Equation  (31) 
used  to  calculate  the  tangent  modulus  drops  drastically  at  high  strain  values,  even  below 
zero  at  especially  high  strains.  Since  a  negative  elastic  modulus  has  no  physical  meaning 
in  the  context  of  the  proposed  model,  this  problem  was  avoided  by  setting  a  lower  limit 
on  the  modulus  of  elasticity.  A  minimum  modulus  of  1 0  MPa,  ~  1 00  times  smaller  than 
the  epoxy  under  zero  strain,  was  chosen  to  capture  the  “failure”  of  the  reinforcement 
without  allowing  catastrophic  failure  of  that  element  and  the  model  as  a  whole. 
Furthermore,  the  present  model  also  accounts  for  the  hollow  structure  of  carbon 
nanotubes.  As  explained  in  the  description  of  the  Embedded  Fiber  Method,  the 
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properties  of  the  reinforcing  fiber  are  directly  combined  with  those  of  the  matrix  material. 
Implicit  in  this  method  is  the  assumption  that  the  entire  volume  of  the  reinforcing  fiber 
contributes  to  the  reinforcement.  In  other  words,  the  reinforcing  fiber  is  assumed  to  have 
a  solid  structure.  First  discovered  by  Iijima  in  1994,  carbon  nanotubes  are  a  novel  form 
of  carbon  and  are  in  fact  hollow  tubes  [68].  The  hollow,  tubular  structure  is  clearly 
shown  in  the  work  of  Pipes  et  al.  [55]  and  is  shown  in  Figure  14. 


Figure  14.  A  diagram  of  a  single  SWCNT  structure  illustrates  the  hollow,  tubular  nature  of  CNTs 
[55]. 

The  assumption,  then,  of  solid  reinforcing  fibers  made  by  the  EFM  is  clearly  not  accurate 
on  the  nanoscale.  This  assumption  would  overestimate  the  reinforcing  ability  of  the 
nanotubes,  if  left  uncorrected,  since  the  stiffness  is  calculated  using  the  cross-sectional 
area  of  the  fiber  as  shown  in  Equation  (21).  Several  authors  have  reported  various 
methods  for  making  this  correction.  Seidel  and  Lagoudas  [41]  developed  a 
micromechanics  model  of  CNT-composites  in  which  they  used  a  composite  cylinders 
approach  to  calculate  effective  properties  of  nanotubes  when  used  as  reinforcing  fibers. 
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Using  an  initial  elastic  modulus  of  1.1  TPa  as  supported  by  the  literature  cited  in  the 
introduction,  they  calculated  an  effective  axial  elastic  modulus  of  704  GPa  by  comparing 
the  strain  energies  of  sold  and  hollow  cylinder  models  of  the  nanotubes.  Odegard  et  al. 
[40]  argue  that  a  standard  micromechanics  approach  is  not  accurate  due  to  the  fact  that 
the  lattice  structures  of  the  nanotube  and  polymer  chains  cannot  be  considered  continuous 
at  the  nanoscale.  They  instead  propose  modeling  the  nanotube/polymer  interface  as  an 
effective  fiber  by  using  an  equivalent-continuum  modeling  method.  They  report  an 
effective  axial  elastic  modulus  of  450.4  GPa.  Considering  this  approach  along  with  the 
small  degradation  of  properties  caused  by  functionalization  [23],  the  fibers  were  assumed 
to  have  elastic  moduli  of  400  GPa  in  the  current  model. 

In  addition  to  the  preceding  data,  a  standard  Poisson’s  ratio  of  w  =  ©3  was  assumed  for 
the  epoxy  matrix.  Nanotube  upper  and  lower  length  limits  of  20  and  800 
respectively,  were  chosen  for  the  angle  of  deviation  calculation  in  Equation  (6).  This 
assumption  is  based  upon  the  reported  length  distributions  of  Wang  et  al.  [29]  and  Ziegler 


et  al.  [30], 
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Chapter  4 

Numerical  Results 

Having  discussed  the  problem  specific  details  and  assumptions  of  the  proposed  approach, 
the  elasticity  and  thermal  results  are  next  presented.  First,  the  deformation  behavior  of 
the  composite  in  tension  is  shown  for  two  different  boundary  conditions.  A  comparison 
of  the  elastic  results  to  the  experiments  detailed  is  then  attempted.  Finally,  a  discussion 
of  the  thermal  results  compared  to  other  models  and  experiments  is  undertaken. 

4.1  Deformation  Behavior 

The  first  result  derived  by  the  current  model  is  a  qualitative  description  of  the 
deformation  behavior  of  the  composite  structure  as  a  whole  for  a  given  applied 
displacement  in  the  positive  y  direction.  By  visualizing  the  displacements  of  the  various 
elements  throughout  the  RVE,  it  is  possible  to  glean  information  about  the  distortion  of 
the  structure  and  the  localized  reinforcement  effects  of  the  carbon  nanotubes.  To 
accomplish  this,  the  characteristics  of  the  epoxy  used  by  Sun  et  al.  [46]  and  described  in 
Chapter  3  were  used.  As  also  indicated  above,  this  is  a  parameter  that  can  be  easily 
altered  based  on  the  particular  composite  matrix  of  interest.  A  reinforcing  component  of 
1  wt%  functionalized  SWCNTs  (F-SWCNTS)  was  randomly  generated  throughout  the 


49 


RVE.  To  conduct  the  deformation  analysis,  the  composite  generated  by  the  model  was 
analyzed  under  two  different  loading  conditions.  In  both  tests,  the  upper  boundary  of  the 
RVE  was  displaced  0.035  pm,  equivalent  of  a  total  effective  strain  of  3.5%  on  the  RVE 
in  the  positive  y-direction. 

In  the  first  test,  the  lower  boundary  of  the  RVE  was  held  fixed  in  the  indirection  only. 
This  test  condition  was  chosen  because  it  eliminates  any  reactions  at  the  fixed  nodes  in 
the  x-direction.  By  doing  so,  the  reactions  at  the  fixed  nodes  can  then  be  used  to  recover 
the  effective  stresses  that  will  be  used  for  the  stress-strain  curve.  A  diagram  of  the 
boundary  conditions  for  the  first  test  condition  is  shown  in  Figure  15. 


Figure  15.  A  diagram  of  the  boundary  conditions  for  the  first  test  condition  is  presented.  The  bottom 
left  corner  is  pinned  and  the  bottom  boundary  is  on  rollers. 

Using  these  boundary  conditions,  the  final  displacements  of  all  the  elements  in  the  finite 
element  mesh  were  examined.  The  number  of  elements  in  the  mesh  is  referred  to  as  the 
“mesh  refinement.”  The  mesh  refinement  necessary  for  experimentally  accurate  results  is 
discussed  in  the  next  section.  A  non-linear  analysis  was  performed  for  one,  ten,  and  one 
hundred  incremental  displacement  steps.  The  different  increment  sizes  were  examined 
for  a  variety  of  different  mesh  refinements.  As  already  mentioned  in  the  description  of 
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the  Embedded  Fiber  Method,  one  of  its  primary  advantages  is  the  capacity  to  derive 
accurate  results  with  a  simple  square  mesh.  Because  the  mesh  is  square,  its  refinement, 
as  reported  in  this  thesis,  will  be  described  using  the  term  “divisions”.  A  division  refers 
to  the  square  root  of  the  number  of  elements  in  the  mesh.  For  example  a  mesh  with  one 
division  has  one  element,  two  divisions  yield  four  elements,  four  divisions  yield  sixteen 
elements,  etc.  This  concept  is  illustrated  in  Figure  16. 


Figure  16.  A  visualization  of  the  mesh  refinement  for  one,  two,  three,  and  five  divisions,  respectively 
is  shown. 

The  final  displacements  of  the  mesh  can  be  visualized  in  two  different  ways.  The  first  is 
to  simply  redraw  the  final  mesh  superimposed  upon  the  initial  mesh.  This  kind  of 
visualization  allows  a  qualitative  view  of  the  shape  change  of  the  RYE.  To  readily 
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observe  this  behavior,  the  incremental  solution  technique  was  applied  with  100  equal 
steps  to  an  RVE  with  10  divisions  and  1  wt%  volume  fraction  of  SWCNTs.  The 
effective  strain  of  the  composite  was  taken  to  3.5%  strain.  The  final  and  initial  meshes 
are  both  shown  in  Figure  17. 

x  10'6 


Figure  17.  The  mesh  deformation  visualization  for  a  single  RVE,  with  10  divisions  and  1  wt% 
SWCNTS,  stretched  to  3.5%  strain  is  presented.  The  change  in  shape  of  the  RVE  can  be  easily  seen 
and  by  examining  the  shape  change  of  each  element,  the  reinforcing  effect  of  fibers  can  be  observed. 
The  displacements  of  the  final  mesh  are  multiplied  by  a  factor  of  three  to  make  visualizing  easier. 


The  second  option  of  visualizing  the  displacements  of  the  mesh  involves  viewing  a 
contour  plot  of  displacements  in  the  positive  y  direction.  This  technique  may  be  more 
visually  suggestive  about  nanotube  reinforcement  as  the  locations  of  the  nanotubes  can  be 
conveniently  shown  along  with  the  contour  plot.  A  contour  plot  of  the  same  RVE  as 
depicted  in  Figure  17  is  now  shown  in  Figure  18. 
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Figure  18.  A  contour  plot  of  the  y  displacements  of  the  element  nodes  is  presented  with  the  nanotube 
locations  superimposed.  The  reinforcing  properties  of  the  nanotubes  can  be  observed  by  noting  the 
regions  where  colors  stay  the  same.  Some  samples  of  these  regions  are  highlighted  with  white  boxes. 

Figure  18  clearly  reveals  the  reinforcing  capabilities  of  the  nanotubes.  The  most  notable 
features  to  observe  in  the  contour  plot  are  the  areas  in  which  color  remains  the  same  or 
nearly  the  same  in  the  y  directions.  These  regions  indicate  quite  small  displacements  and 
are  clearly  located  where  nanotubes  are  aligned  in  the  tensile  stretch  direction.  A  few  of 
these  regions  are  outlined  in  white  for  easy  identification. 

Each  of  these  techniques  to  visualize  displacement  is  useful  alone,  but  when  they  are 
viewed  together,  a  comparison  can  be  readily  made,  enhancing  the  descriptive  ability  of 
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both  plots.  By  viewing  the  two  plots  side  by  side,  the  reinforcing  effects  of  the  nanotubes 
can  be  seen  not  only  in  the  tensile  direction,  but  also  in  the  x  direction.  To  this  effect,  the 
mesh  deformation  and  displacement  contours  of  the  same  RVE  are  shown  in  Figure  19. 
However,  in  this  figure  the  plots  are  shown  for  a  single  step  at  a  high  strain  value, 
approximately  2.5%  strain.  By  examining  the  deformation  behavior  at  high  strains,  the 
reinforcing  capabilities  are  more  readily  observable  since  the  surrounding  matrix  material 
is  losing  its  capacity  to  elastically  resist  the  tensile  stretching. 


Figure  19.  The  mesh  deformation  and  contour  plot  of  y  displacements  for  the  same  RVE  are 
presented  at  2.5%  strain.  By  visualizing  the  displacements  of  a  single  incremental  step  the 
reinforcing  capabilities  of  the  nanotubes  are  very  clear.  The  displacements  of  the  final  mesh  on  the 
left  are  multiplied  by  300  to  make  visualizing  easier. 


Furthermore,  it  is  insightful  to  compare  the  contour  plots  of  the  RVE  before  and  after  the 
onset  of  plastic  deformation.  By  making  this  comparison,  it  becomes  evident  that  at  low 
strains,  the  nanotubes  have  a  limited  reinforcing  role,  but  at  high  strains  their  role 
becomes  very  important  and  pronounced.  The  contour  plots  of  the  displacements  after 
the  first  step  and  after  the  75"'  step  are  shown  in  Figure  20. 
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Figure  20.  The  contour  plots  of  the  y  displacements  after  the  first  and  75th  incremental  steps  are 
presented.  By  examining  the  deformation  before  and  after  the  onset  of  plastic  deformation,  the  role 
of  the  nanotubes  in  combating  plastic  deformation  is  evident. 


Having  presented  the  deformation  results  from  the  first  loading  condition  as  depicted  in 
Figure  15,  the  results  of  the  second  test  condition  are  now  shown.  The  second  condition 
is  the  traditional  tensile  test  in  which  the  bottom  edge  of  the  sample  is  fixed  in  both  the  x 
and  the  y  directions  and  the  top  edge  is  fixed  in  the  x  direction  while  displaced  in  the  y 
direction.  These  boundary  conditions  are  shown  in  Figure  21. 


Figure  21.  A  diagram  of  the  boundary  conditions  for  the  second  test  condition  illustrates  the 
traditional  tensile  test.  The  bottom  edge  is  completely  fixed  and  the  top  edge  is  fixed  in  the  x- 
direction. 
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In  the  same  manner  as  the  first  test  condition,  the  change  in  the  RVE  shape  can  be 
visualized  by  the  change  in  the  mesh  shape.  The  final  and  initial  shapes  are  shown  for  an 
RVE  with  no  nanotubes  and  an  RVE  with  1  wt%  SWCNTs  for  comparison  in  Figure  22. 
In  these  results,  100  incremental  steps  were  again  used  but  20  divisions  were  used  this 
time  as  compared  to  the  10  divisions  used  in  the  results  above.  This  difference  is  made 
only  to  show  what  20  divisions  look  like  when  visualized. 


Figure  22.  The  mesh  deformation  for  the  second  testing  condition  is  presented,  (a)  Mesh  deformation 
for  an  RVE  without  nanotubes.  Symmetric  deformation  is  observed  (b)  Mesh  deformation  for  an 
RVE  with  1  wt.  %  SWCNTs  is  shown.  The  final  mesh  displacements  are  multiplied  by  a  factor  of  3 
to  make  visualizing  easier. 
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Similarly,  the  contour  plots  with  and  without  nanotubes  are  shown  in  Figure  23. 


Figure  23.  The  displacement  contours  for  the  second  test  condition  are  presented,  (a)  The 
displacement  contour  of  RVE  without  reinforcement  is  shown,  (b)  The  displacement  contour  of  an 
RYE  with  1  wt  %  nanotubes  is  shown  for  comparison. 


An  enlargement  of  the  contour  plot  in  Figure  23  (b)  is  shown  in  Figure  24.  Here,  the 
areas  outlined  in  white  highlight  the  reinforcing  capability  of  the  nanotubes.  It  is  also 
worth  recognizing  the  path  in  which  major  changes  in  color  occur.  These  paths,  indicated 
by  white  lines,  clearly  follow  regions  that  are  void  of  nanotubes  indicating  a  lack  of 


mechanical  reinforcement. 
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Figure  24.  A  displacement  contour  of  the  second  test  condition  with  special  features  indicated  is 
presented.  The  areas  boxed  in  white  indicate  sample  regions  of  particularly  strong  reinforcement  in 
the  tensile  direction  as  shown  by  the  continuity  of  color.  The  white  paths  across  the  RVE  indicate 
probable  locations  of  fracture  in  the  RVE  due  to  lack  of  reinforcement. 

After  using  both  testing  conditions  to  examine  deformation  behavior,  the  first  condition 
was  chosen  to  calculate  the  stress-strain  curve  for  the  composite.  As  detailed  above,  by 
allowing  all  the  nodes  to  move  freely  in  the  x  direction  as  the  first  condition  does,  there 
are  no  reaction  forces  in  the  x  direction.  Hence,  all  of  the  force  caused  by  displacing  the 
upper  boundary  nodes  will  be  present  in  the  y  component  of  the  reaction  forces  at  the 
bottom  nodes.  This  fact  allows  a  true  stress-strain  relationship  for  the  lone  composite 
material  to  be  developed. 
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4.2  Stress-Strain  Results 

With  a  qualitative  view  of  the  deformation  behavior  in  hand,  the  next  step  was  to 
determine  the  mesh  refinement  and  number  of  incremental  steps  needed  for  an  accurate 
model.  To  accomplish  this,  a  mesh  refinement  analysis  was  performed  and  the  results 
were  compared  to  experimental  data.  The  stress-strain  curve  used  for  comparison  in  the 
mesh  refinement  analysis  is  the  functionalized-SWCNT  curve  (top,  blue  curve)  of  Sun  et 
al.  [46]  shown  in  Figure  1 1 . 

Due  to  the  random  nature  of  the  nanotube  geometry  generation,  the  finite  element  results 
of  a  single  RVE  would  not  be  a  statistically  accurate  representation  of  the  stress-strain 
relationship  of  the  general  composite.  To  report  a  statistically  accurate  result,  a  Monte 
Carlo  statistical  analysis  must  be  adopted  for  the  data.  Monte  Carlo  statistical  analysis  is 
a  process  in  which  meaningful  information  about  a  random  event,  such  as  the  stress-strain 
curve  of  any  single  RVE  in  the  current  model,  is  gained  by  statistically  analyzing  a  large 
number  of  those  events.  In  the  case  of  the  current  model,  the  only  analysis  that  is  needed 
is  to  take  the  average  stress  values  of  many  different  RVEs  at  the  same  strain  values  to 
create  a  statistically  accurate  stress-strain  curve.  To  ensure  that  an  adequate  number  of 
random  RVE’s  were  considered,  the  non-linear  simulation  is  performed  on  500  different 
RVE  microstructures.  The  average  stress  values  for  each  step  are  then  used  to  generate 
the  stress-strain  curve  for  the  composite  material.  Esteva  [42]  validated  the  choice  of  five 
hundred  different  microstructures  by  performing  a  Monte  Carlo  convergence  analysis  on 
his  linear  model.  Since  the  proposed  non-linear  model  is  a  combination  of  many  linear 
solutions,  this  analysis  validates  five  hundred  samples  as  a  sufficient  sample  size  for  the 


current  model  as  well. 
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Before  reporting  the  stress-strain  curve  for  a  material,  it  was  first  necessary  to  determine 
the  number  of  incremental  steps  and  the  number  of  divisions  required  for  the  non-linear 
solution  to  converge  to  an  accurate  value.  Using  the  same  volume  fraction,  nanotube  and 
epoxy  densities,  and  Poisson’s  ratio  as  Sun  et  al.  [46],  the  Monte  Carlo  analysis  described 
above  was  performed  on  the  results.  To  determine  the  number  of  incremental  steps 
needed,  the  results  of  the  incremental  solution  for  a  single  division  of  the  RVE  were 
examined.  This  convergence  analysis,  along  with  the  F-SWCNT  curve  of  Sun  et  al.  [46] 
for  reference,  is  summarized  in  Figure  25. 
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Figure  25.  (a)  The  non-linear  incremental  convergence  analysis  for  the  model  is  shown  for  a  single 
division  of  the  RVE.  (b)  A  closer  view  of  the  model  results  at  high  strains.  A  step  refinement  of  100 
incremental  steps  is  identified  as  sufficient  for  convergence. 
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By  examining  these  results,  100  incremental  steps  were  identified  as  necessary  to  capture 
the  non-linearity  of  the  epoxy.  The  Monte  Carlo  analysis  was  then  performed  on  the 
results  of  500  different  RVE’s  for  1,10,  20,  40,  60,  and  80  divisions  in  order  to  determine 
the  mesh  refinement  needed  for  an  accurate  result.  This  kind  of  mesh  convergence 
analysis  is  herein  termed  Monte  Carlo  Mesh  Convergence  Analysis  (MCMMA).  The 
MCMMA  for  100  incremental  steps  is  presented  in  Figure  26  along  with  the  F-SWCNT 
curve  reported  by  Sun  et  al.  [46]. 
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Figure  26.  (a)  The  mesh  convergence  analysis  for  the  model  with  one  hundred  incremental  steps  is 
shown  for  1,  10,  20,  40,  60,  and  80  divisions  of  the  RVE.  The  “diamond”  line  identifies  the 
experimental  results  of  Sun  et  al.  [45].  (b)  A  closer  view  of  the  model  results  at  high  strains.  A  mesh 
refinement  of  60  divisions  is  identified  as  sufficient  for  accurate  stress-strain  results. 
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The  mesh  convergence  analysis  shown  above  indicates  that  an  accurate  result  can  be 
obtained  with  60  divisions  of  the  finite  element  mesh.  Though  a  finer  mesh  will  improve 
the  accuracy,  as  it  can  be  observed  by  the  line  for  80  divisions,  the  computational  cost 
increases  as  n 2,  where  n  is  the  number  of  divisions.  This  fact  creates  the  need  to  identify 
the  minimum  number  of  divisions  necessary  for  accurate  results,  and  is  shown  to  be  60 
divisions  in  Figure  26. 

To  confirm  the  choice  of  100  incremental  steps  and  make  sure  that  accurate  results 
cannot  be  obtained  by  simply  increasing  the  number  of  divisions,  the  same  MCMMA  as 
shown  in  Figure  26  was  conducted  for  both  one  and  ten  incremental  steps.  These  results 
are  shown  in  Figure  27.  It  is  readily  seen  that  neither  one  nor  ten  steps  yield  accurate 
results  for  any  number  of  divisions. 
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Figure  27.  (a)  The  mesh  convergence  analysis  for  the  mode!  with  one  incremental  step  is  shown  for  1, 
10,  20,  40,  60,  and  80  divisions  of  the  RVE.  The  “diamond”  line  identifies  the  experimental  results  of 
Sun  et  al.  [45].  (b)  The  same  analysis  is  shown  for  the  model  with  10  incremental  steps. 
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With  the  values  of  60  divisions  and  100  incremental  steps  selected  as  the  model 
parameters  needed  for  qualitatively  accurate  results,  the  final  step  was  to  compare  the 
results  of  the  model  to  different  experimental  results  in  order  to  validate  the  accuracy  and 
effectiveness  of  the  model.  For  these  data  sets,  statistical  analysis  was  performed  to 
determine  the  standard  deviation  of  the  Monte  Carlo  data.  The  standard  deviation  of  a 
given  data  set  is  determined  by  the  equation  [69] 

a  =  -  x)2  .  (34) 

In  this  equation,  N  is  the  number  of  data  points,  x,  is  an  individual  data  point,  and  x  is  the 
mean  value  of  that  data  set.  After  calculating  the  standard  deviation  of  the  stresses  at  the 
different  strain  levels,  the  stress-strain  curve  for  60  divisions,  100  incremental  steps  was 
compared  to  two  different  experimental  data  sets 

The  stress-strain  results  of  the  model  were  first  compared  to  the  experimental  results  of 
Sun  et  al.  [46],  The  tangent  modulus  of  the  epoxy  matrix  was  updated  using  Equation 
(30)  and  the  tensile  test  was  simulated  up  to  3.5%  strain.  The  comparison  is  shown  in 
Figure  28. 
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Figure  28.  (a)  The  model  results  for  60  divisions,  100  steps  are  compared  to  the  results  of  Sun  et  al. 
[45].  (b)  A  close  view  of  the  stress-strain  curve  high  strain  values  is  shown. 
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As  can  be  observed  Figure  28,  the  experimental  curve  is  close,  but  does  not  fall  within 
the  standard  deviation  error  bars  of  the  model  at  all  levels  of  strain.  The  maximum 
disagreement  between  the  model  and  the  experiment  is  less  than  4%.  Because  the 
experimental  curve  is  derived  from  a  single  realization  of  an  experiment,  this  lack  of 
agreement  is  not  an  indication  of  a  poor  model.  By  looking  at  the  curves  of  the  minimum 
and  maximum  RVE’s  shown  in  Figure  29,  a  large  span  of  values  for  the  stress  at  high 
strain  levels  can  be  observed.  The  experimental  curve  of  Sun  et  al.  falls  within  this 
region.  This  fact  indicates  that  the  model  is  in  fact  a  good  one  and  the  lower  curve 
reported  by  Sun  et  al.  could  be  accounted  for  by  any  number  of  factors. 


60  Div,  100  Steps,  Min  and  Max  RVE's 


Figure  29.  The  maximum  and  minimum  RVE’s  out  of  the  500  RVE’s  sampled  in  the  MCMCA  are 
compared  to  the  experimental  results  of  Sun  et  al.  [45]  to  highlight  the  span  of  possible  outcomes  for 
random  nanotube  geometries. 


In  addition  to  the  results  of  Sun  et  al.,  the  model  is  also  compared  to  the  experimental 
data  reported  by  Zhu  et  al.  [22],  The  equation  used  to  update  the  elastic  modulus  for  this 
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case  is  detailed  in  Equation  (33).  In  this  case,  the  model  is  used  for  up  to  5%  strain  and 
once  again,  60  divisions  and  1 00  incremental  steps  were  used.  The  comparison  of  the 
model  to  Zhu  et  al.’s  data  is  shown  in  Figure  30. 


Figure  30.  The  model  results  for  60  divisions,  100  steps  are  compared  to  the  results  of  Zhu  et  al.  [22], 
The  center,  blue  F-SWCNT  line  represents  the  experimental  data. 


Akin  to  the  first  experiment,  the  model’s  calculated  stress-strain  curve  does  not  fall 
within  the  standard  deviation  error  bars  at  all  strain  values.  By  simply  increasing  the 
mesh  refinement  of  the  model,  a  closer  curve  could  easily  be  generated.  However,  the 
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goal  of  the  model  proposed  herein  is  to  merely  yield  accurate  results,  and  this  is 
accomplished  with  60  divisions  as  demonstrated  above.  The  largest  disagreement 
between  the  model  and  the  experimental  data  is  less  than  5%.  Because  of  the  close 
agreement  with  two  different  experiments,  the  elastic  model  is  deemed  reliable.  Utilizing 
a  non-linear  finite  element  procedure  with  only  60  divisions  in  a  simple  square  mesh,  an 
accurate  stress-strain  curve  can  be  generated  for  an  arbitrary  nanotube-reinforced 
polymer  composite  if  the  polymer  matrix’s  stress-strain  curve  is  known. 

4.3  Thermal  Results 

The  final  result  of  the  approach  presented  in  this  thesis  is  calculation  of  the  effective  a 
non-linear  thermal  conductivity  of  the  composite.  Both  the  incremental  and  iterative 
non-linear  approaches  were  considered  in  the  thermal  model  to  determine  the  superiority 
of  one  over  the  other  in  terms  of  computational  time  and  numerical  accuracy.  For  the 
thermal  case,  only  the  non-linear  thermal  conductivity  of  SWCNTs  was  considered 
because  its  scale  is  at  least  four  orders  of  magnitude  greater  than  the  thermal  conductivity 
of  the  polymer  matrix,  at  any  temperature.  The  thermal  conductivity  of  the  polymer  used 
for  the  present  model  is  assumed  to  be  constant  because  of  the  experimental  evidence  of 
both  Biercuk  et  al.  [50]  and  Xu  et  al.  [47],  First,  the  same  method  used  to  visualize  the 
displacement  contours  in  the  mechanical  analysis  is  applied  to  the  thermal  problem  to 
visualize  the  temperature  contour  of  any  RVE.  Mesh  refinement  analysis  was  then 
performed  on  the  results  of  both  the  incremental  and  iterative  approaches  and  the 
effective  thermal  conductivity  for  130  divisions  was  compared  to  the  values  reported  in 
the  current  literature.  For  the  mesh  refinement  analysis,  a  polymer  thermal  conductivity 
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of  0.188 - -  was  used.  This  value  is  in  close  agreement  with  the  values  reported  by 

Wang  et  al.  [48]  for  their  epoxy  resin,  and  Cai  and  Song  [49]  for  their  polyurethane 
matrix. 

Boundary  Conditions  and  Temperature  Contours 

To  calculate  the  effective  thermal  conductivity  of  the  composite,  the  non-linear  thermal 
properties  of  SWCNTs  were  taken  from  the  work  of  Grujicic  et  al.  [15],  A  graph  of  the 
non-linear  dependence  of  thermal  conductivity  on  temperature  is  presented  in  Figure  9. 
In  this  model,  the  bottom  edge  of  the  RVE  is  held  at  0°  C  ,  the  top  edge  is  held  at  100°  C, 
and  both  side  edges  are  insulated.  As  usual,  these  values  could  be  altered  depending  on 
the  experimentalist’s  interests.  The  thermal  conductivity  of  the  composite  was  calculated 
using  the  heat  flux  at  the  bottom  nodes  of  the  RVE.  A  diagram  of  these  boundary 
conditions  is  shown  in  Figure  3 1 . 
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T=  100  °C 

v  v 


Figure  31.  A  diagram  of  the  boundary  eonditions  for  the  thermal  analysis  is  shown.  The  top  edge  is 
held  at  100°  C  and  the  bottom  is  held  at  0°C  while  the  sides  are  insulated. 


To  compare  to  the  elasticity  analysis,  a  given  displacement  was  applied  to  the  top  of  the 
RYE  and  the  reaction  forces  at  the  bottom  fixed  nodes  were  recovered  and  used  to 
calculate  the  stress  in  the  RVE.  For  the  thermal  analysis,  a  temperature  difference  was 
applied  between  the  top  and  bottom  of  the  RVE  and  the  heat  fluxes  at  the  bottom  nodes 
were  recovered.  The  thermal  conductivity  of  the  RVE  was  then  calculated  using  the 
equation 


K  — 

K«ff  ~ 


(35) 


where,  E  q  is  the  sum  of  the  heat  fluxes  across  the  bottom  nodes,  A  is  the  area  of  the 
RVE,  t  is  the  thickness  of  the  RVE,  and  AT  is  the  temperature  difference  between  the 


bottom  and  top  of  the  RVE. 
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Before  calculating  the  thermal  conductivity,  the  first  result  of  the  thermal  analysis  was  a 
qualitative  view  of  the  thermal  conductivity  enhancement  provided  by  the  addition  of 
SWCNTs  to  the  polymer.  With  a  similar  motivation  to  the  deformation  behavior  section 
of  the  mechanical  analysis,  visualizing  a  temperature  contour  of  the  RVE  allows  the  user 
to  gain  a  good  understanding  of  how  effectively  nanotubes  conduct  heat  through  the 
composite.  Using  the  polymer  thermal  conductivity  of  0.188  -  and  the  nanotube 

thermal  conductivity  curve  of  Grujicic  et  al.  [15],  the  non-linear  thermal  model  was  run 
using  10  iterative  steps  and  40  divisions.  As  will  be  seen  in  the  following  section,  the 
number  of  steps  is  more  than  enough  to  capture  the  non-linearity  of  the  nanotubes  and  the 
number  of  divisions  is  sufficient  to  simply  see  the  temperature  contours  accurately.  The 
temperature  contour  for  a  random  RVE  with  1  wt%  SWCNTs  is  shown  in  Figure  32. 
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Figure  32.  A  contour  plot  of  temperatures  from  the  non-linear  thermal  model  is  presented. 

Similar  to  the  displacement  contours  for  the  elastic  model,  the  temperature  contour 
vividly  illustrates  the  reinforcing  effect  of  the  nanotubes.  The  areas  of  the  contour  plot 
where  color  remains  constant  indicate  regions  of  high  thermal  conductivity.  These  areas 
clearly  follow  the  nanotubes  in  the  RVE  and  indicate  their  role  in  enhancing  the  thermal 
conductivity  of  the  composite. 
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Mesh  Convergence  Analysis 

After  achieving  a  qualitative  confirmation  that  the  model  was  working  correctly  through 
the  temperature  contours,  the  next  step  was  to  determine  the  number  of  steps  necessary 
for  convergence  in  the  non-linear  solution  techniques.  For  both  the  iterative  and 
incremental  approaches,  the  thermal  conductivities  of  the  nanotubes  were  updated  after 
each  step  using  the  temperatures  of  the  nanotubes  in  conjunction  with  Equation  (29).  The 
final  thermal  conductivity  of  the  composite  is  calculated  using  the  heat  flux  at  the  last 
step  in  conjunction  with  Equation  (35).  Monte  Carlo  Mesh  Convergence  Analysis 
(MCMCA)  was  used  to  obtain  statistically  accurate  values  of  thermal  conductivity  and  to 
determine  the  finite  element  mesh  size  necessary  for  the  thermal  conductivity  to  converge 
to  an  accurate  value.  500  RVE’s  were  considered  again  for  the  same  reasons  outlined  in 
the  stress-strain  results  section. 

The  first  approach  taken  to  consider  the  non-linear  thermal  conductivity  of  the  nanotubes 
was  the  incremental  approach  described  in  Chapter  3.  Temperature  increments  between 
0°  C  and  100°  C  were  applied  to  the  top  boundary  of  the  RVE.  The  MCMCA  for  the 
incremental  approach  is  shown  in  Figure  33. 
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Figure  33.  (a)  The  MCMCA  for  the  thermal  conductivity  of  500  RVEs  is  presented  for  the 
incremental  approach.  A  mesh  refinement  of  130  divisions  is  needed  for  accurate  results,  (b)  By 
zooming  in,  it  is  also  evident  that  10  incremental  steps  are  sufficient  to  capture  the  non-linearity  of 
the  thermal  conductivity  of  the  nanotubes. 
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The  MCMCA  shown  in  Figure  33  highlights  the  three  important  results  of  the  thermal 
model.  First,  an  effective  thermal  conductivity  of  ~  0.80  — —  was  calculated  for  the 

composite  with  1  wt.  %  SWCNT  loading.  This  result  shows  that  the  model  is  indeed 
yielding  accurate  thermal  conductivities.  Second,  a  mesh  refinement  of  130  divisions 
was  identified  as  necessary  for  the  FEA  to  converge  to  a  single  value.  Finally,  it  was 
determined  that  10  incremental  steps  are  sufficient  to  capture  the  non-linear  thermal 
conductivity  of  the  SWCNTs. 

The  second  approach  taken  to  capture  the  non-linearity  of  the  SWCNTs  was  an  iterative 
approach.  In  this  approach,  a  single  temperature  gradient  of  1 00°  C  was  applied  to  the 
RVE.  Using  the  resultant  temperatures  throughout  the  RVE,  the  thermal  conductivities 
of  the  nanotubes  were  updated  and  the  gradient  is  applied  again.  MCMCA  was  used  to 
investigate  the  number  of  steps  needed  to  capture  the  non-linearity  of  the  nanotubes  as 
well  as  the  mesh  refinement  needed  for  a  convergent  effective  thermal  conductivity  of  the 
composite.  This  analysis  is  summarized  in  Figure  34. 
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Figure  34.  (a)  The  MCMCA  for  the  thermal  conductivity  of  500  RVEs  is  presented  for  the  iterative 
approach.  A  mesh  refinement  of  130  divisions  is  needed  for  accurate  results,  (b)  By  zooming  in,  it  is 
also  evident  that  only  two  iterative  steps  are  needed  to  capture  the  non-linearity  of  the  thermal 
conductivity  of  the  nanotubes. 
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The  MCMCA  shown  in  Figure  34  highlights  the  three  important  results  of  the  thermal 
model  as  well.  First,  an  effective  thermal  conductivity  of  ~  0.80  —  was  calculated  for 

the  composite  with  1  wt.  %  SWCNT  loading.  Second,  a  mesh  refinement  of  130 
divisions  was  identified  as  necessary  for  the  FEA  to  converge  to  a  single  value.  Finally, 
it  was  determined  that  2  iterative  steps  are  sufficient  to  capture  the  non-linear  thermal 
conductivity  of  the  SWCNTs. 

Comparing  these  two  different  approaches,  it  is  apparent  that  for  a  mesh  refinement  of 
130  divisions,  both  the  incremental  and  iterative  approaches  yield  almost  identical 
results.  However,  it  takes  the  incremental  approach  10  steps  to  accomplish  the  same  non¬ 
linear  analysis  that  takes  only  two  steps  in  the  iterative  approach.  For  this  reason,  the 
iterative  approach  is  chosen  as  the  preferred  non-linear  approach  to  take  for  the  present 
model. 

Comparison  to  Experiment 

The  final  and  most  important  step  in  the  thermal  analysis  of  composites  using  the  current 
model  is  the  comparison  of  the  results  to  experimental  data.  After  a  thorough  review  of 
the  recent  literature  on  the  subject,  the  data  of  Wang  et  al.  [48],  Moisala  et  al.  [16], 
Biercuk  et  al.  [50],  Hong  and  Tai  [51],  Cai  and  Song  [49],  and  Xu  et  al.  [47]  was  chosen 
to  provide  a  broad  variety  of  polymers  and  volume  fractions  for  comparison.  Agreement 
with  experimental  results  both  validated  the  thermal  model  and  confirmed  its  utility. 
Before  these  comparisons  are  made,  it  should  be  first  noted  that  it  was  often  difficult  to 
compare  experiments  to  the  present  model  exactly  for  several  reasons.  First, 
experimental  reports  do  not  always  report  the  complete  properties  of  the  composite 
matrix  that  they  used.  The  most  common  property  left  out  is  the  density  of  the  matrix. 
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This  prevents  a  truly  accurate  conversion  of  weight  fraction  to  volume  fraction.  An 
approximate  density  of  1.2  ^7  was  used  in  these  cases  to  convert  weight  fractions  to 

volume  fractions  using  Equation  (7).  Second,  different  kinds  of  nanotubes  and 
procedures  are  often  used.  For  example,  Xu  et  al.  [47]  reported  results  for  non-purified 
SWCNTs  with  no  functionalization  while  Wang  et  al.  [48]  used  shortened  nanotubes  of 
50  rmi  or  less  that  are  well  dispersed  throughout  the  matrix.  Moisala  et  al.  [16]  and  Cai 
and  Song  [49]  reported  results  for  MWCNT’s.  Hong  and  Tai  [51]  used  thin  films  of 
SWCNT  networks,  called  Buckypapers,  in  reporting  their  results.  In  short,  there  are 
various  experimental  approaches  taken  to  report  thermal  conductivity  of  these 
composites,  and  while  the  model  proposed  herein  can  account  for  a  wide  variety  of 
physical  parameters  of  the  matrix  and  nanotubes  themselves,  it  cannot  account  for  all  of 
the  intricacies  mentioned  above.  However,  as  the  MCMCA  results  above  showed,  one 
can  expect  the  model  to  still  yield  reliable  values  for  the  thermal  conductivities  of 
composites. 

First,  the  model  was  compared  to  six  different  experiments.  For  each  comparison,  the 
model  was  executed  using  the  same  volume  fraction  and  matrix  thermal  conductivity, 
as  reported  in  the  papers  from  which  the  data  was  taken.  These  comparisons  are 


shown  in  Figure  35. 
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Figure  35.  A  comparison  of  the  thermal  conductivities  calculated  by  the  model  to  various 
experiments  is  shown. 


As  is  obvious  from  Figure  33,  the  proposed  model  results  are  close  to  relevant 
experimental  results  for  small  weight  fractions  of  1  wt%  or  less.  However,  it  is  also 
evident  that  the  model  overestimates  the  thermal  conductivity  of  the  composite  in  every 
single  case.  Both  of  these  observations  can  be  explained  by  the  imperfect  assumption  of 
perfect  transfer  of  properties  between  the  nanotubes  and  the  polymer  matrix  inherent  in 
the  Embedded  Fiber  Method. 

The  primary  reason  for  the  inaccuracy  of  this  assumption  relates  to  the  high  thermal 
resistance  between  the  nanotube-matrix  and  nanotube-nanotube  interfaces.  Huxtable  et 
al.  [70]  found  that  the  interface  conductance,  G,  for  CNT’s  in  a  polymer  matrix  is 
approximately  12  .  For  comparison,  Wilson  et  al.  [71]  found  an  interface 


81 


conductance  of  130  -h—  for  platinum  nanoparticles  in  water.  This  exceptionally  small 

interface  conductance  will  result  in  a  significant  reduction  in  the  thermal  conductivity  of 
the  composite.  Closely  related,  Xue  [72]  performed  a  study  examining  the  effects  of 
interfacial  resistance,  or  Kapitza  resistance,  between  the  nanotubes  and  the  surrounding 
matrix  on  the  effective  conductivity  of  nanotube-reinforced  composites.  His  results 
indicate  that  this  resistance  can  effectively  reduce  the  thermal  conductivity  of  the 
nanotubes  by  several  orders  of  magnitude.  While  the  small  thermal  conductance  and 
high  thermal  resistance  accounts  for  the  small  overestimation  of  thermal  conductivity  at 
low  volume  fractions,  it  does  not  adequately  explain  the  discrepancy  observed  at  higher 
volume  fractions. 

The  significant  disagreement  between  the  model  and  experiments  for  high  volume 
fractions  can  be  explained  by  the  large  number  of  junctions  among  carbon  nanotubes  that 
exist  as  the  nanotube  volume  fraction  increases.  While  this  percolating  network  is  ideal 
for  electrical  conductivity,  a  large  resistance  to  phonon  movement  from  one  nanotube  to 
another  exists  for  thermal  loads  [47],  This  resistance  actually  dominates  the  small 
interfacial  conductance  and  is  the  primary  reason  why  dramatic  increases  in  thermal 
conductivity  are  seen  for  small  volume  fractions  compared  to  very  small  gains  with  the 
addition  of  more  nanotubes.  For  example,  Wang  et  al.  [48]  report  a  40%  increase  in 
thermal  conductivity  of  their  composite  with  the  addition  of  only  0.5  wt%  SWCNTs.  In 
contrast,  Xu  et  al.  [47]  report  only  a  60%  increase  in  thermal  conductivity  with  19% 
volume  fraction  of  SWCNTs. 

With  rather  limited  research  reported  in  these  two  areas  of  thermal  conductivity 
hindrance,  it  is  impossible  to  incorporate  these  factors  empirically  into  the  model. 
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However,  the  phenomenological  observation  of  the  composite  thermal  conductivities 
reported  in  the  literature  has  allowed  the  incorporation  of  these  factors  into  the  model  by 
simply  reducing  the  thermal  conductivity  of  the  SWCNTs  by  a  factor  of  twenty.  While 
still  using  the  same  non-linear  thermal  conductivity  reported  by  Grujicic  et  al.  [15],  the 
result  of  Equation  (29)  is  simply  divided  by  twenty  and  then  the  rest  of  the  model 
continues  normally.  Using  this  observation,  the  effective  thermal  conductivities  shown  in 
Figure  32  were  recalculated  and  again  compared  to  the  same  experiments.  This 
comparison  is  shown  in  Figure  36  which  reflects  a  much  better  agreement  between  the 
model  results  and  the  experimental  data. 
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Figure  36.  A  comparison  of  the  modified  model’s  thermal  conductivity  results  to  various  experiments 
is  presented.  Much  better  agreement  with  experimental  results  is  clearly  evident. 

The  comparison  shown  in  Figure  34  clearly  illustrates  that  the  simple  reduction  of  the 

thermal  conductivity  of  the  nanotubes  by  a  factor  of  twenty  creates  much  better 

agreement  between  the  model  and  experiment  for  all  volume  fractions.  For  every 
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comparison,  except  that  to  Xu  et  al.,  the  thermal  conductivity  predicted  by  the  model  is 
within  ~  33%  of  the  value  reported  by  the  corresponding  experiment.  While  this  value 
may  seem  large  at  first,  considering  the  poorly  understood  and  complex  mechanisms 
described  above  that  inhibit  efficient  heat  conduction,  this  agreement  is  indeed  quite 
reasonable. 

In  conclusion,  the  thermal  conductivity  calculated  using  a  two-step  iterative  finite 
element  procedure  using  the  Embedded  Fiber  Method  produced  results  that  provide  a 
good  approximation  of  the  thermal  properties  of  a  SWCNT-reinforced  polymer 
composite  for  volume  fractions  less  than  about  three  percent.  Comparing  model  results 
to  six  different  related  experiments,  and  finding  good  agreement,  demonstrates  the 
versatility  of  the  model  for  simulating  a  variety  of  different  polymers  over  a  large  range 
of  SWCNT  volume  fractions.  However,  the  accuracy  of  the  model  could  be  improved 
with  a  better  understanding  of  the  mechanisms  inhibiting  efficient  heat  transfer  between 
the  polymer  and  the  nanotubes,  and  between  nanotube  junctions. 
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Chapter  5 


Concluding  Remarks 


In  the  past  decade,  a  significant  amount  of  experimental  and  theoretical  research  has  been 
conducted  to  develop  and  characterize  a  large  number  of  novel  nanocomposites.  The 
experimental  research  has  been  significantly  limited  by  the  cost,  complexity,  and  time 
required  in  developing  and  testing  new  nanocomposites.  To  date,  no  theoretical  model 
has  been  able  to  incorporate  the  non-linear  mechanical  and  thermal  properties  of  these 
new  composites  for  a  sample  of  significant  reinforcing  volume  fraction. 

In  this  thesis,  a  non-linear  finite  element  model  has  been  presented  to  address  this  need. 
For  mechanical  properties,  the  model  was  utilized  to  visualize  the  deformation  behavior 
and  to  calculate  the  stress-strain  curve  of  a  SWCNT-reinforced  polymer  composite  under 
tensile  strain.  For  thermal  properties,  the  model  presented  was  used  to  calculate  the 
effective  thermal  conductivity  of  the  composite  by  incorporating  the  non-linear  thermal 
conductivity  of  the  SWCNTs. 

The  three  major  advantages  of  the  proposed  model  have  been  exemplified  in  the  results 
presented  in  this  paper.  First,  the  capacity  to  visualize  the  deformation  of  a  SWCNT- 
reinforced  polymer  was  demonstrated  through  plots  of  the  original  and  deformed  finite 
element  meshes  along  with  contour  plots  of  the  nodal  displacements  throughout  the  RYE. 
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By  comparing  these  two  visualizations,  along  with  the  positions  of  nanotubes  in  the  RVE, 
significant  information  about  the  reinforcing  effects  of  the  nanotubes  can  be  gleaned.  In 
addition,  the  temperature  contour  plot  presented  in  the  thermal  results  allows  the 
enhancement  in  thermal  conductivity  of  the  composite  due  to  the  nanotubes  to  be  clearly 
observed. 

Second,  the  ability  to  generate  an  accurate  stress-strain  curve  of  the  composite  has  been 
demonstrated.  Using  a  mesh  refinement  of  60  divisions  along  with  a  100  incremental 
steps,  the  non-linear  elastic  properties  of  the  polymer  matrix  were  accurately  modeled. 
After  comparing  the  model  results  to  the  work  of  Sun  et  al.  [46]  and  Zhu  et  al.  [22],  the 
model  has  been  found  to  be  accurate  within  5%  for  the  entire  stress-strain  curve. 

Finally,  the  capacity  to  calculate  the  effective  thermal  conductivity  of  a  SWCNT- 
reinforced  polymer  composite  was  shown.  By  incorporating  the  non-linear  thermal 
properties  of  SWCNTs,  the  accuracy  of  the  model  is  improved.  Using  a  mesh  refinement 
of  130  divisions  along  with  2  iterative  steps,  the  thermal  conductivities  of  seven  different 
experimental  scenarios  have  been  calculated.  After  comparing  the  results  of  the  model  to 
the  work  of  Wang  et  al.  [48],  Moisala  et  al.  [16],  Biercuk  et  al.  [50],  Hong  and  Tai  [51], 
Cai  and  Song  [49],  and  Xu  et  al.  [47],  it  has  been  found  that  the  model  is  accurate  within 
33%  for  any  volume  fraction  less  that  3%.  While  this  error  may  seem  large,  it  is  in  fact 
quite  reasonable  when  considering  the  poorly  understand  mechanisms  inhibiting  heat 
transfer  between  CNTs  and  between  CNTs  and  the  polymer. 

In  summary,  it  has  been  found  that  a  non-linear  finite  element  model  using  the  Embedded 
Fiber  Method  can  accurately  predict  the  mechanical  and  thermal  properties  of  a  large 
variety  of  different  SWCNT-reinforced  polymer  composites.  This  model  has  many 
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potential  uses  in  composite  materials  engineering.  It  provides  experimentalists  with  a 
tool  to  validate  their  results,  as  well  as  a  means  to  accurately  explore  approximate 
properties  of  systems  for  future  experimentation.  The  primary  feature  of  this  model  is 
its  flexibility  for  use  in  a  variety  of  applications.  The  model  proposed  in  this  thesis 
provides  a  timely  and  accurate  calculation  of  deformation,  stress-strain  relationships,  and 
thermal  conductivities  for  any  number  of  new  composite  materials. 

Further  work  in  the  direction  of  this  thesis  may  include  the  calculation  of  electrical 
properties  for  these  SWCNT-reinforced  polymer  composites,  modifications  to  allow  the 
modeling  of  DWCNT  and  MWCNT-reinforced  polymer  composites,  and  the  extension  of 
the  model  to  three  spatial  dimensions. 
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Appendix  A 


Newton  Rhapson  Scheme  to  Find  the  Natural  Coordinates  of  a  Nanotube 

subroutine  N  atCoordF  inder(exc,x,y  ,nc) 

!  find  natural  coordinates  from 
!  element  global  coordinates  (reverse  mapping) 
implicit  none 

real(8),  intent(in)  ::  exc(8,l),x,y 
real(8),intent(out)  ::  nc(2,l) 
real(8)  xl,x2,x3,x4,yl,y2,y3,y4 
real(8)  al,a2,a3,a4,bl,b2,b3,b4 
real(8)  J(2,2),Ji(2,2),f(2,l),nci(2,l),djj 

!  extracting  node  coordinates 

xl=exc(l,l) 

yl=exc(2,l) 

x2=exc(3,l) 

y2=exc(4,l) 

x3=exc(5,l) 

y3=exc(6,l) 

x4=exc(7,l) 

y4=exc(8,l) 

al  =0.25*(xl  +x2+x3+x4) 
a2=0.25*(-xl-x2+x3+x4) 
a3=0.25*(-xl+x2+x3-x4) 
a4=0.25*(xl-x2+x3-x4) 

b  1  =0 .2  5  *  (y  1  +y2+y  3 +y4) 
b2=0 .2  5  *  (-y  1  -y2+y  3  +y4) 
b3=0.25*(-yl+y2+y3-y4) 
b4=0 .2  5  *  (y  1  -y2+y  3  -y4) 

nc=0.d0  !  initial  guess 
d=l.d0 

do  while  (d>0.0001d0) 

!  Calculating  Jacobian  matrix 
J(  1 , 1  )=a3+a4*nc(2, 1 ) 

J(l,2)=a2+a4*nc(l,l) 

J(2, 1  )=b3+b4*nc(2, 1 ) 

J(2,2)=b2+b4*nc(  1,1) 

!  Calculating  vector  of  equations 
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f(  1 , 1  )=a  1  +a2  *nc(2, 1  )+a3  *nc(  1 , 1  )+a4*nc(  1 , 1 )  *nc(2, 1  )-x 
f(2, 1  )=b  1  +b2*nc(2, 1  )+b3  *nc(  1 , 1  )+b4*nc(l ,  1  )*nc(2, 1  )-y 
!  Calculating  inverse  of  Jacobian  matrix 
jj=J(l,l)*J(2,2)-J(l,2)*J(2,l) 

Ji(l,l)=J(2,2)/jj 
Ji(l,2)=-J(l,2)/jj 
Ji(2, 1  )=-J(2,  l)/jj 
Ji(2,2)=J(l,l)/jj 

!  Newton-Rhapson  Approximation 
nci=nc-matmul(Ji,f) 

!  Calculating  the  norm  of  difference 
d=sqrt((nci(  1 , 1  )-nc(  1,1))**  2+(nci(2, 1  )-nc(2, 1 ))  *  *2) 
nc=nci 
enddo 


end  subroutine  NatCoordFinder 


